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DIGEST 


The  pa£>c  decade  has  witnessed  the  rapid  assembly  of  rather  large 
body  of  data  relating  to  the  mechanical  behavior  and  failure  of  materials 
under  extreme  pressure,  temperature  and  shear  loading  conditions.  Al¬ 
though  this  has  led  to  a  number  of  important  phenomenological  insights, 
the  molecular  engineering  of  new  materials  with  improved  response  and 
strength  characteristics  in  impeded  by  the  absence  of  a  fundamental  under¬ 
standing  of  the  failure  modes  of  substances  subjected  to  these  conditions. 

For  the  past  several  years  a  group  at  the  Vitreous  State  Laboratory 
of  Catholic  University  of  America  has  been  engaged  in  the  investigation  of 
the  response  of  liquids  and  amorphous  solids  to  large,  rapidly  imposed 
disturbances.  Initially  our  approach  was  to  make  use  of  Molecular  Dynamics 
(MD)  experiment  simulations  to  probe  the  microscopic  origins  and  mechanisms 
of  certain  anomalous  macroscopic  response  features  of  liquids  under  such 
conditions. 

In  its  initial  stages  the  work  was  focused  on  such  non-linear  visco¬ 
elastic  phenomena  as  shear  viscosity  thinning  and  stress  overshoot  under 
conditions  of  large  shearing  rates.  In  carrying  out  this  work  it  became 
apparent  that  the  behavior  that  was  observed  in  the  simple  MD  model  system 
(of  point  particles  interacting  via  a  Lennard-Jones  6-12  potential)  were 
quite  general — having  been  seen  in  such  diverse  systems  as  organic  lubri¬ 
cants,  polymers,  metals,  glasses  and  composites.  In  each  of  these  cases, 
the  extant  explanations  of  the  behavior  involved  mechanisms  specific  to  the 
particular  structural,  configurational,  or  bonding  properties  of  the  indi¬ 
vidual  systems;  the  MD  data  suggest  that  the  origin  of  many  of  these  obser¬ 
vations  can  be  understood  in  a  more  general  fundamental  way. 

In  particular,  our  studies  showed  that  it  is  possible  to  rationalize 
much  of  the  behavior  as  a  consequence  of  the  strain-induced  (or  stress-in¬ 
duced)  structural  changes  that  were  observed  to  occur  in  the  MD  experiments. 

In  the  liquid  and  high  temperature  glass  systems,  for  example,  it  was  found 
that  large  shear  rates  cause  the  development  of  locally  organized  anisotro¬ 
pic  structures  in  terms  of  which  the  dynamical  pseudo-plastic  behavior  of 
these  materials  can  be  understood.  These  structural  changes  result  in  a 
limit  to  the  steady-state  stress  that  the  material  can  support  and  offer  a 
framework  for  understanding  the  ductile  failure  process  in  these  substances. 

To  complement  the  MD  work  a  supporting  experimental  program  was  estab¬ 
lished  with  in-house  university  funds.  This  joint  experiment al/MD  effort 
has  proven  to  be  unusually  valuable  both  in  assisting  in  the  interpretation 
of  the  data  that  is  obtained,  and  in  providing  guidance  for  continuing  studies 
It  was  found,  for  instance  that  effective  viscosity  data  from  both  the  MD 
and  laboratory  experiments  exhibited  essentially  the  same  dependence  on  im- 
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posed  serein  rate  and  could  be  interpreted  in  terms  of  the  same  non¬ 
linear  viscoelastic  response  hypothesis.  The  analysis  of  the  MD  results 
showed  that  the  crucial  material  property  in  the  data  reduction  can  be 
represented  by  a  parameter  that  can  be  interpreted  as  the  maximum  steady- 
state  stress  that  can  be  developed  in  the  material  under  the  application 
of  a  uniform  strain  note.  In  the  experimental  work  it  was  found  that  when 
a  material  is  subjected  to  a  stress  exceeding  this  value,  it  initially  de¬ 
forms  plastically;  however,  as  the  stress  is  continued,  this  ductile  re¬ 
sponse  is  arrested  and  the  material  fails  in  sudden  fracture. 

This  transition  of  ductile  flow  to  fracture  occurs  as  the  material 
attempts  to  rearrange  its  structure  to  reduce  the  stress  accompanying  the 
applied  strain  rate.  When  these  structural  changes  cannot  occur  rapidly 
enough  to  compensate  for  the  development  of  stress  levels  in  excess  of  the 
limiting  stress  (initially  conceived  as  he  parameter  used  in  fitting  the 
MD  and  experimental  data),  the  material  fractures  spontaneously.  The  lim¬ 
iting  steady-state  stress  value  Xs  therefore  the  actual  cohesive  strength 
of  the  material. 

In  the  MD  studies  of  low-temperature  glasses  the  imposition  of  large 
strain  rates  causes  initially  a  pseudo-plastic  material  response  leading  to 
ductile  failure;  under  confinement  conditions  the  material  then  restructure 
itself  in  a  crystalline  form  with  improved  modulus  and  strength  properties. 
Some  preliminary  studies  in  a  crystalline  system  suggest  that  qualitatively 
similar  behavior  (the  development  of  transient  dislocations  followed  by  re¬ 
crystallization)  occurs  at  elevated  shear  stress  levels  in  crystals. 


INTRODUCTION 


This  technical  report  presents  the  results  of  research  in 
the  dynamical  response  of  materials  to  large,  suddenly  applied 
disturbances.  The  work  reported  herein  was  carried  out  under  the 
sponsorship  of  the  Office  of  Naval  Research  Research  (Contract 
Numbers  N00014-75-C-08S6  and  N00014-81-K-0296)  and  covers  the 
time  period  October  1,  1979  through  May  31,  1981. 

The  body  of  the  report  consists  of  six  sections.  In  the 
first  is  presented  a  molecular  dynamics  (MD)  investigation  'f  the 
non-linear  shear  response  of  a  Lennard-Jones  model  system  at 
elevated  pressures.  Among  the  key  results  reported  here  are  the 
observations  of  "shear  thinning"  and  "shear  overshoot"  effects  at 
high  shear  rates.  These  are  accompanied  by,  and  indeed  are,  a 
consequence  of  dynamical  strain-induced  organization  of  the 
liquid  structure.  In  the  second  section,  these  data  are  com¬ 
pared  with  results  obtained  in  both  inorganiz  glass  and  polymeric 
systems.  While  there  are  quantitative  differences  among  the 
various  materials,  several  important  features  are  qualitatively 
quite  similar,  indicating  that  the  molecular  dynamics  studies 
(where  one  is  able  to  extract  detailed  microscopic  information) 
can  be  quite  useful  in  the  interpretation  of  phenomena  observed 
in  "real  materials."  The  third  section  continues  this  comparison 
of  the  molecular  dynamics  "experimental  results"  with  primary 
laboratory  data.  In  this  work  it  is  shown  that  the  concept  of 
shear  induced  structural  organization  derived  from  the  MD  work  can 
be  applied  in  characterizing  and  understanding  the  mechanical  and 
thermal  responses  of  organic  traction  lubricants  at  high  shearing 
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rates . 

To  interpret  these  results  a  model  proposed  by  Bair  and 
Winer  [S.  Bair  and  W.O.  Winer,  ASME  J.  Lub.  Tech.  100 ,  40, 

(1979)]  was  generalized  to  provide  a  mathematical  framework 
in  wnich  to  consider  dynamical  nonlinear  viscoelasticity.  This 
constitutes  the  fourth  section  of  the  report.  One  of  the  key 
features  in  this  work  is  the  identification  of  the  dimensionless 
quantity 

x  *  ^t0(Gco/°*) 

(here  i  =  strain  rate,  t  =  shear  relaxation  time,  G  =  shear 

O  oo 

rigidity  modulus,  and  a*  =  maximum  steady-state  stress)  as  the 
parameter  governing  the  type  of  behavior,  i.e.,  linear  visco¬ 
elastic,  pseudoplastic,  etc.,  that  is  observed.  It  was  hypothe¬ 
sized  that  for  large  x  values  ductile  failure  of  the  material 
will  be  the  result.  In  sections  five  and  six,  MD  and  primary 
laboratory  results  obtained  in  such  circumstances  are  given, 
respectively.  In  these  experiments  low  temperature  glassy 
materials  for  which  tq  (and  therefore  x)  is  large  were  studied 
and  the  predicted  failure  was  indeed  observed.  In  the  MD  work 
the  initial  pseudoplastic  response  of  the  material  gave  way  to 
ductile  failure  in  which  the  stress  that  was  developed  dropped 
precipitously  to  zero.  Under  confinement  conditions  the  continued 
application  of  shear  led  to  strain-induced  crystallization  as  the 
material  restructured  itself  in  a  higher  strength  form. 


Section  I. 


Time  dependent  non-linear  shear  stress  effects  in  simple 
liquids:  A  molecular  dynamics  study. 


by  D.  M.  He  ye  3 ,  J 
and  T. 


J«  Kim,  C.  1 
A.  Litovitz 


.  Mont  rose 


Tima  dependent  nonlinear  shear  stress  effects  In  simple 
liquids:  A  molecular  dynamics  study** 

O.  M.  Heyes,  J.  J.  Kim,bl  C.  J.  Montrose,  and  T.  A.  Litovitz 

furrow  Start  laboratory.  Catholic  Unvarsity  of  America,  Washington  D.C.  100M 
(Received  IT  December  1979;  accepted  1 1  July  1910) 

TV  effect*  of  larye  amplitude  shearing  rate*  on  104  particle  amorphous  Lcnnard-Jones  systems  have  been 
examined.  A  resolution  of  the  structural  and  dynamical  properties  reveals  that  there  is  a  tendency  for  shear  to 
reorganite  the  liquid  into  layers  to  facilitate  flow  and  hence  reduce  shear  viscosity.  The  first  steps  have  been 
made  to  determine  those  physical  phenomena  that  are  necessary  tn  develop  a  theory  of  time  dependent 
nonlinear  shear  viscoelastic  effects  The  form  of  the  distortion  giving  rise  to  nonlinear  shear  stresses  is  well 
understood  although  its  magnitude  is  intimately  linked  with  parallel  structural  changes  in  an,  at  present, 
intractable  way. 


I.  INTRODUCTION 

In  a  rather  wide  variety  o(  technological  applications, 
e.g. ,  elastohydrodynamic  lubrication,  shock  loading, 
fatigue,  liquids  and  amorphous  solids  are  subjected  to 
large  shearing  forces.  The  microscopic  mechanisms 
that  are  of  importance  in  such  nonlinear  situations  are 
not  well  understood.  To  attempt  to  gain  some  insights 
into  this  phenomenon,  a  molecular  dynamics  (MD)  in¬ 
vestigation  of  the  viscoelastic  respo.tse  of  liquids  sub¬ 
jected  to  large  shear  strain  rates  was  undertaken. 

We  have  used  a  modification  of  the  normal  MD  pro¬ 
cedure  in  order  to  follow  the  response  of  the  system  to 
an  applied  perturbation.  In  the  past,  nonequilibrium 
MD  has  been  used  to  produce  steady  shear  flow1  in  a 
number  of  ways.  As  hurst  and  Hoover2  sheared  a  molec¬ 
ular  system  through  the  action  of  “fluid  walls"  on  two 
opposite  faces  of  the  MD  cell.  Their  technique  produces 
surface  effects  which  lead  to  an  inhomogeneous  material 
being  sheared.  This  problem  is  largely  eliminated  by 
the  homogeneous  shear  method3  (HSM),  which  uses  pe¬ 
riodic  boundary  conditions  in  all  directions.  A  modifi¬ 
cation  of  the  HSM,  similar  to  one  proposed  by  Evans, 4 
was  used  in  these  calculations.  It  enables  a  particular 
strain  rate  to  be  established  instantaneously.  The  un¬ 
perturbed  velocity  of  the  real  or  image  molecule  i,  i.e. , 
v(,  is  altered  to  Vt  by  applying  a  shear  rate  e„  so  that 
it  only  alters  the  x  component  of  velocity 


where  is  the  ^  coordinate  of  particle  i  measured  from 
the  center  of  the  cubic  MD  cell  of  side  length  L.  A 
molecule  that  leaves  the  MD  cell  in  the  t  direction  is  re¬ 
introduced  through  the  opposite  face  with  the  (x,.v)  coor¬ 
dinates  and  velocity  of  its  nearest  displaced  imap*.  All 
the  calculations  were  conducted  isothermally  by  scaling 
the  v,  to  maintain  the  desired  temperature. 


*’ Research  supported  in  part  by  the  office  of  Naval  Reaearoh 
Contract  No.  NOOOU-75-C-0856. 

“On  leave  from  Korea  Advanoed  Institute  of  Science,  p.  o. 
Box  ISO  Chonguangni,  Seoul,  Korea. 


Calculations  were  performed  under  steady  state  shear 
conditions.  In  addition,  the  response  of  the  system  to  a 
time  dependent  strain  rate  was  studied.  These  latter 
calculations  consisted  of  segments  lasting  for  80S  time 
steps.  In  each  segment  the  same  “large”  strain  rate 
c„  was  applied  during  the  first  lime  step  and  was  main¬ 
tained  for  the  next  430  steps,  at  which  time  it  was  sud¬ 
denly  switched  off  and  then  held  at  xero  for  the  remain¬ 
der  of  the  segment.  The  time  dependent  system  prop¬ 
erties  were  averaged  over  approximately  20  such  seg¬ 
ments  in  order  to  reduce  the  statistical  noise. 

The  segmented  MD  procedure  was  modified  to  allow 
the  response  of  the  system  to  small  shear  strain  rates 
5e„  to  be  followed. 5  The  system  was  allowed  to  evolve 
from  the  same  starting  configuration  with  and  without 
the  perturbation.  A  great  deal  of  the  statistical  noise 
was  eliminated  by  subtracting  at  each  time  step  the 
value  for  the  property  under  examination  of  the  undis¬ 
turbed  system  from  that  of  the  perturbed  medium. 

The  MD  cell  contained  108  particles  interacting 
through  the  Lonnard-Jones  (LJ)  potential 

<p(r)  =  4«[(«r/r)12 -(o/r)8] .  (2) 

The  units  of  energy,  length,  and  time  used  throughout 
this  paper  are  c ,  o,  and  (mo*/<  )*’,  respectively;  m  is 
the  particle’s  mass.  All  other  quantities  are  expressed 
in  terms  of  these  fundamental  units.  The  equations  of 
motion  were  integrated  using  the  Verlet  algorithm8  with 
a  time  step  having  a  typical  duration  of  0.005.  The  in¬ 
teractions  were  truncated  beyond  »■=  2. 5o.  The  system 
was  investigated  at  an  average  temperature  T  of  0. 722 
(c/*j)  and  number  densities  p  of  0.8442,  0.92862, 
1.01304,  and  1.2663.  Had  each  of  the  latter  three  states 
been  achieved  slowly,  the  system  would  have  crystal¬ 
lized7;  however,  these  states  were  in  fact  achieved  by 
suddenly  (in  one  time  step)  reducing  the  volume  of  the 
p  =  0.8442  system  by  10%,  30%,  and  50%,  respectively. 
This  rapid  denslficatlon  prohibited  those  structural  re¬ 
arrangements  that  are  necessary  for  crystallization  to 
occur  so  that  a  "superdenalfied”  liquid  system  was  pro¬ 
duced.  Accumulation  of  the  configurational  averages 
followed  an  initial  equilibration  period  of  several  thou¬ 
sand  time  steps.  The  calculations  were  performed  on  a 
PDP  11/40  computer. 
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Heyet  et  § I . :  Nonlinear  shear  in  liquids 
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FIG.  1.  The  time  dependence  of  the  reduced  sheer  atreaa 
cu(»)/C.  for  the  p«  1.01304  systems  subjected  to  the  follow ln( 
strain  rate  history:  cw(f)«c,  0<H2.3.  <„(!!■  0,  2.3<f(4.6; 
•••*,  £*0.1793; - ,  £  *0.3584;  end  — ,  4*0.7168. 


II.  RESULTS  AND  DISCUSSION 

We  first  discuss  the  nature  of  the  relaxational  pro¬ 
cesses  that  accompany  a  time  varying  large  strain  rate. 
The  segmented  technique  vaa  used.  Each  segment  con¬ 
sisted  of  continuous  shear  for  0  « t  2. 3.  The  system 
was  then  allowed  to  evolve  in  the  absence  of  shear  from 
1*2.3  to  1*4.  5,  when  the  shear  rate  was  again  applied. 

The  development  of  the  stress  tensor  o  Is  of  particular 
interest;  the  afi  component  of  this  was  computed  using 
the  formula* 


(3) 

where  V  is  the  volume  of  the  MD  cell  which  contains  N 
molecules;  i>*  is  the  a  component  of  the  velocity  of  the 
1th  perticle  relative  to  the  Imposed  velocity;  atl  and  plt 
are  the  o  and  0  components,  respectively,  of  the  vector 
r(Jar,  -  tj,  where  r,  defines  the  position  of  molecule  i. 

Immediately  after  a  atraln  rate  « ,,  is  applied  at  time 
1  =  0,  a  liquid  takes  time  to  respond  structurally.  As  a 
result,  at  first  the  shear  stress  rises  linearly  with  time 
as  in  an  elastic  material 

o,*(f)  =  G,  cr,f  .  (4) 

where  Gm  is  the  shear  rigidity  modulus.  Since  a  liquid 
will  not  support  a  static  stress  (it  will  flow),  the  rate  of 
ascent  of  shear  stress  decreases  until  it  reaches  a 
limiting  value  governed  by  the  shear  vlscorlty  r>,  i.e. , 
o„(l  -  •)  = »)  £  „.  This  is  illustrated  in  Fig.  1,  which 
shows  the  shear  rate  dependence  of  o„(f)/G„  for  the  p 
=  1 . 01304  state.  Observe  that  the  steady  state  stress 
does  not  increase  in  proportion  to  the  shear  rate,  i.e., 
the  shear  viscosity  decreases  with  increasing  shear 
rate— so-called  shear  thinning.  Table  I  reveals  that 
this  trend  is  manifest  at  all  densities.  This  is  strong 
evidence  that,  on  shearing,  structural  changes  are  in¬ 
duced  in  a  liquid  which  alters  its  physical  properties. 

A  second  experiment  was  carried  out  to  probe  the 
viscosity  behavior  at  high  shear  rates.  For  a  given  c„, 


TABLE  I.  The  density  p  and  ateady  state  shear  rate  dependence  of  the  shear  visooaity  t), 
normal  preaeura  components  oonflgurational  energy  U„  shear  rigidity  modulus  G„,  and 
atreaa-optloal  coefficient  c(.  Values  of  G,  obtained  from  Eq.  (3)  are  given  in  parentheses. 
The  estimated  maximum  standard  error  la  *  10%  at  the  lowest  shear  rate  of  each  density. 


p 

e 

n 

Pm 

Pr, 

Pm 

v. 

G. 

0.8448 

0.0 

—  2.70 

0.8 

0.8 

0.8 

-6.88 

22.5 

(23.7) 

0.8448 

0.0422 

2.64 

0.8 

0.9 

0.9 

-5.64 

0.8443 

0.1686 

2.58 

1.0 

0.9 

1.0 

-5.62 

0. 8442 

0.3373 

2.44 

1.2 

1.0 

1,3 

-5.58 

0. 8442 

0.8745 

2.12 

1.6 

1.4 

1.7 

-5.61 

3.5 

0. 92882 

0.0 

-10 

3.0 

3.0 

3.0 

-6.05 

33,4 

(37.3) 

0. 92882 

0.0435 

7.4 

3.0 

2.9 

3.0 

-6.05 

0.92862 

0. 1741 

4.6 

3.3 

3.2 

3.4 

-5.99 

0. 92862 

0.3482 

4.0 

3.7 

3.5 

5.9 

-5.91 

0. 92862 

0.8963 

3.3 

4.4 

4.2 

4.7 

-5.77 

4.0 

1.01304 

0.0 

-20 

<7.0 

<7.0 

<7.0 

<-6.38 

46.9 

(49.4) 

1.01304 

0.0898 

12.6 

7.0 

8.4 

7.0 

-8.22 

4.0 

1.01304 

0. 1792 

9.0 

7.5 

7,4 

7.3 

-6.14 

4.0 

1.01304 

0. 3584 

5.7 

8.3 

7.9 

8.4 

-5.96 

4.5 

1.01304 

0.7188 

4.6 

9.  S 

8.8 

9.9 

-5.75 

6.0 

1.2863 

0.0 

-200 

<37 

<37 

<37 

<-4.6 

123.8 

(108.6) 

1.2683 

0.0483 

104 

40 

36 

38 

-4.49 

1.2883 

0.1930 

25 

42 

39 

41 

-3.94 

1.2883 

0.3881 

16 

48 

43 

48 

-3.29 

1.2683 

0.7721 

11 

49 

47 

48 

-2.68 

8.6 
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I  Hi,  2.  The  time  dependent  normal  pressure  components  lor 
the  strain  rate  history  of  Fig.  i  and  t  -  0.  710H,  p  -  I.  0130-4. 

— .  - i  snd  Pmtt). 


the  system  wan  allowed  to  reach  its  steady  state  stress 
value  c„,  at  which  time,  say  t',  a  small  increment  in 
shear  rate  fic„  was  applied.  The  stress  8 ct,(t)  result¬ 
ing  from  the  application  of  8<u  continuously  from  t'  is 
given  by  linear  response  theory  as 

8o„(f)  =  G-  b(  x,(t")C(t -t")  ,  (5) 

•V 

where  the  relaxation  function  to  be  used  in  Eq.  (5)  is 
that  characteristic  of  the  sheared  state  (see  below). 

This  leads  directly  to  an  alternate  measurement  of  the 
shear  viscosity  at  high  shear  rates: 

T)  =  lim8c„(<)/8<,(  .  (6) 

»-• 

The  values  obtained  in  this  manner,  although  less  pre¬ 
cise,  agree  reasonably  well  *  1th  and  follow  qualitatively 
the  same  trends  as  those  show),  in  Table  I  derived  from 

The  time  evolution  of  the  shear  a  ress  is  also  affected 
rather  dramatically  at  the  higher  shear  rates.  In  Fig. 

1  it  can  be  seen  that  the  shear  strain  rises  to  a  peak 
before  descending  to  its  steady-state  value.  This  be¬ 
havior  is  observed  for  shear  rates  such  that  cjj/C. 

>  0. 05;  at  lower  shear  rates  no  local  maxima  in  the  o>( 
versus  t  curves  were  discernible.  One  also  finds  that 
upon  cessation  of  the  applied  shear  rate,  the  immediate 
stress  relaxation  is  faster  for  the  more  highly  sheared 
media.  The  rate  of  stress  decay  slows  with  time  as  the 
equilibrium  state  is  approached.  We  shall  later  consider 
the  possible  microscopic  origin  of  the  shear  thinning  and 
shear  stress  overshoot  phenomena;  before  doing  so, 
however,  we  examine  some  of  the  other  changes  that  oc  - 
cur  in  the  system  in  response  to  the  large  and  suddenly 
applied  shear  disturbances. 

Each  normal  pressure  component  Pum  -  -  omm  rises 
when  the  liquid  is  sheared.  They  reflect  the  structural 
changes  that  the  shear  forces  introduce  Into  the  liquid. 
The  P«„(f)  of  the  p=  1.01304,  c«, =0. 7168  state  are  given 
in  Fig.  2.  The  time  scale  of  the  development  of  the  nor¬ 


mal  pressure  components  Is  longer  than  that  of  the  shear 
stress.  No  normal  pressure  overshoot  is  present  and  a 
P„  overshoot  is  observed  only  when  1.2063.  Table 
1  reveals  that  the  configurational  energy  per  particle 

H  V 

Ut  r,  (i  /2N)Y.  T.  <?>(>•,,) 

I*/ 

increases  with  shear  rate.  This  is  additional  evidence 
of  structural  alteration  of  the  fluid  to  states  which  are 
onh  energetically  favorable  in  the  presence  of  steady 
shear  flow. 

The  shear  stress  viscoelastlcic  behavior  la  usefully 
cast  In  relaxation  function  formaliam.  The  normalised 
shear  stress  relaxation  function  C(f)  has  been  evaluated 
from  both  the  sheared  and  unsheared  media  using  the 
aforementioned  small  perturbation  technique.  The 
change  in  shear  stress  6o„(/)  following  the  application 
of  a  small  strain  rate  6c„  for  the  duration  of  one  time 
step  At  at  /'  is  given  by  C(t  - /')  as  follows: 

8c„(f  -  tf)  =  Gmt>iIIAtC(t  - 1')  ,  (7) 

where  C(Q)  =  1  and  C’(*>)  =  0.  The  C(t  - 1')  derived  for  the 
p=  1. 01304  state  with  c„=  0. 7168,  taking  time  origins 
at/'=:0. 16,  1.2,  and  3.6,  are  shown  in  Fig.  3.  At 
short  times  (t  - 1'  sO. 2),  all  the  C(t  - 1 ')  coincide,  pre¬ 
sumably  because  this  portion  of  stress  decay  is  inertial 
and  results  from  small  motions  of  molecules  within  the 
free  apace  of  the  random  network  structure.  After  this, 
the  remaining  portion  of  the  stress  decays  more  slowly 
because  It  Involves  the  structural  rearrangement  of  the 
medium,  requiring  the  participation  of  the  cooperative 
motions  of  larger  molecular  groups.  Figure  3  reveals 
that  as  shear  stress  increases,  the  more  slowly  decay¬ 
ing  process  progressively  disappears.  This  has  been 
observed  experimentally*  and  is  further  evidence  of  the 
structural  changes  that  take  place  In  the  liquid  at  high 
shear  stresses  to  make  shear  flow  easier.  For  the  p 
=  0. 8442  state,  the  form  of  C(<)  with  «„  =  0  Is,  within 
statistical  uncertainty,  the  same  as  that  given  by 
Levesque,  Verlet,  and  KurW'trvl. 10  The  variation  for 
nonzero  values  of  e„  is  qualitatively  the  same  as  for  the 
p  =  1.01304  state,  although  the  changes  in  the  long-time 
behavior  are  somewhat  less  pronounced. 


FIG.  3.  The  shear  e'reae  relaxation  (unction  C(t-t')  for  the 

state  of  Fig.  2  and  time  origin  f  at - ,  f  »0.  t6i  — ,  f 

■  l,2t  and  ••••  l'»3.8. 
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Fill.  I,  The  shear  stress  relaxation  (unctions  C(t)  for  the 

sheared  ( — )  and  unsheared  ( - 1  portions  of  the  segment  for 

the  state  p  - 0.  928(52,  «„•=  0.6963. 


A  selection  of  the  C(/)  obtained  at  the  other  densities 
in  the  steady  state  shear  and  zero  shear  portions  of  the 
segment  are  presented  in  Figs.  4  and  5.  Their  forms 
are  consistent  with  the  above  interpretation.  It  is  cer- 
*ainty  not  surprising  that  the  characteristic  response 
functions  are  different  for  the  highly  stressed  steady 
state  system  as  compared  with  the  equilibrium  system. 
On  fairly  general  grounds,  one  expects  that  the  existence 
of  large  shear  stresses  in  the  liquid  reflects  the  fact 
that  local  structural  distortions  are  present.  That  these 
should  impact  on  the  response  characteristics  of  the 
liquid  seems  quite  apparent. 

The  values  obtained  for  the  shear  rigidity  modulus  C„ 
are  presented  in  Table  I,  where  they  are  compared  with 
the  values  derived  using  the  expression  of  Zwunzig  and 
Mountain11: 

G.  =  pkBT  +  — ^ p2  jT  drifo(r)jp(r4d<p/dr)  .  (8) 

In  this  equation  g0(r)  is  the  radial  distribution  function, 
which  was  computed  as  a  part  of  the  MD  runs  at  each 
density.  Within  the  statistical  uncertainties,  G.  was 
found  to  be  independent  of  the  shear  stress  at  a  particu¬ 
lar  density. 


C(t) 


t 

PIG.  5.  The  shear  stress  relaxation  functions  Of)  lor  the 

sheared  (— )  and  the  unsheared  ( - )  portions  of  the  segment 

for  the  state  p«  1.2663,  f„»0.7T21. 


I-’IG.  6.  The  excess  directional  number  density  functions  6nfa) 
for  the  p-  1.2663  and  *,,  =  0.7721  state. 


We  now  turn  to  consider  the  nature  of  the  structural 
and  dynamiccl  changes  that  these  large  shear  rates  pro  - 
duce  in  the  LJ  liquids.  For  the  purpose  of  investigating 
shear-induced  structural  anisotropy,  a  directional  prob¬ 
ability  density  function  n(a)  was  defined:  pL2n{ot)da 
=  the  probability  of  locating  a  particle  in  the  plane  lamina 
(of  thickness  da  and  dimensions  LxL)  oriented  normal 
to  the  a  axis  and  located  a  perpendicular  distance  a 
from  a  given  particle.  Structural  anisotropy  is  then  re¬ 
vealed  by  examining  6 n(a),  the  difference  between  >’s(ct) 
characterizing  the  sheared  system  and  n0(a)  of  the  equi¬ 
librium  system.  In  terms  of  the  pair  distribution  func¬ 
tion  this  can  be  written  as  (for  the  a  =  z  case,  for  ex¬ 
ample) 

6w(z)  =  (l/L2)  /  dx  f  dy  [&,(r) -^(r)]  ,  (9) 

•/o  yo 

where  g,(r)  is  the  pair  distribution  function  for  the 
sheared  system  and  g0(r)  is  the  equilibrium  radial  dis¬ 
tribution  function.  We  have  obtained  the  6 n(a)  by  com¬ 
puting  the  average  number  of  particles  from  a  reference 
particle  in  the  distance  range  Ixl  to  I x  1  +  Ax  (where  Ax 
was  taken  as  L/200)  selected  from  the  N  - 1  particles  in 
thfc  reorientated  MD  cell.  Figure  6  shows  the  6 n(a)  un¬ 
der  the  conditions  p  =  1. 2663  and  «„=  0. 7721.  Although 
the  maximum  displacement  Is  109c  at  most,  there  is 
evidence  of  a  tendency  of  the  molecules  to  reorder  them¬ 
selves  into  xy  layers  (layers  normal  to  the  velocity  gra¬ 
dient).  The  a'z  shearing  increases  the  probability  of  the 
molecules  in  an  xz  plane  being  found  at  multiples  of  the 
intermolecular  diameter  from  each  other  in  the  z  direc¬ 
tion.  It  is  energetically  favorable  for  these  layers  to  be 
staggered  in  the  y  direction.  Although  not  so  well  de¬ 
fined,  these  trends  are  present  at  the  lower  densities. 
The  (grossly  simplified)  picture  that  emerges  is  one  of 
corrugated  sheets  of  atoms  sliding  past  one  another.  A 
pictorial  representation  of  this  suggested  structure  is 
given  in  Fig.  7. 
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FIG.  8.  ’Itie  trajectories  of  five  molecules  in  the  xz  plane 
during  the  p  =  l. 01304  and  £M  =  0.71G8  segment.  The  molecules 
originated  from  within  n  xz  slice  of  thickness  ~  L/l 0  taken 
through  the  center  of  the  cell. 


0. 


dt 


KIG.  7.  A  pictorial  representation  of  the  structural  changes 
that  take  place  on  going  from  an  unsheared  (a)  to  a  sheared  (b) 
state.  The  arrows  denote  the  line  and  plane  of  shear. 


The  structural  reorganization  results  from  the  action 
of  the  imposed  velocity  profile  which,  at  the  shear  rates 
considered,  dominates  the  molecular  motion.  This  la 
well  demonstrated  by  the  trajectory  plots  of  Fig.  8, 
which  show  time -elapsed  positions  in  the  xz  plane  of  a 
selection  of  the  particles  from  the  pc:  1, 01304  and 
-  0.  7168  state.  Major  movements  are  confined  to  the 
line  of  shear. 

It  is  to  be  expected  that  the  effect  of  the  shear  should 
not  be  confined  to  structural  changes  but  must  also  alter 
the  setf-diffasional  dynamics  of  the  molecules.  To 
assess  this,  the  x,  y,  and  z  components  of  the  mean 
square  displacements  (r„(/)),  which  exclude  the  imposed 
flow,  have  been  calculated,  i,  e. , 

=  E  ([jf  »•?«')<«'  -O|(0)]?)  /N  .  (10) 

The  directional  diffusion  coefficients 


TAULE  11.  The  density  and  shear  rate  depen¬ 
dence  of  the  directional  diffusion  coefficients  Da 
obtained  from  the  MD  experiments  using  Eq.  (10). 
The  results  are  taken  from  the  unsheared  and 
sheared  portions  of  the  segments. 


p 

■D, 

0, 

D. 

0.8442 

0.0 

0.032 

0.032 

0.032 

0.8442 

0.6745 

0.052 

0.045 

0,053 

0.  92862 

0.0 

0.018 

0.018 

0.018 

0.  92862 

0.6963 

0,  043 

0.037 

0.045 

1.01304 

0.0 

0.008 

0.008 

0.008 

1.01304 

0.1792 

0.02 

0.014 

0.02 

1.01304 

0.3684 

0.03 

0.025 

0.03 

1 . 2663 

0.0 

0.0045 

0.0045 

0.0045 

^266^ 

0.7721 

0.045 

0.03 

0.03 

for  a  selection  of  the  states  are  presented  in  Table  II. 
The  value  of  Da  «  0.  032  for  the  unsheared  liquid  near  the 
triple  point  compares  favorably  with  the  value  of  0. 033 
previously  obtained.10  The  D„  of  the  unsheared  liquid 
decrease  by  approximately  an  order  of  magnitude  on  in¬ 
creasing  the  density  from  p-0.8442  to  1.2663.  How¬ 
ever,  at  larger  shear  rates  the  density  has  a  diminish¬ 
ing  influence  on  the  diffusion  coefficients.  The  high 
shear  stresses  enhance  the  fluidity  and  create  a  '‘liquid” 
which  has  directionality.  All  the  diffusion  coefficient 
components  rise  with  shear  rate  at  a  particular  density 
because  the  structural  reorganization  that  takes  place  on 
shearing  apparently  creates  paths  along  which  particles 
can  more  readily  move.  The  are  typically  two  thirds 
of  the  Ox  and  D„  indicating  that  the  self -diffusion  is 
favored  in  the  shearing  plane.  Thus,  paradoxically,  al¬ 
though  structural  aspects  of  the  liquids  under  shear 
present  a  more  solid-like  appearance,  other  more  dy¬ 
namically  related  properties  such  ac  self -diffusion  and 
shear  viscosity  manifest  changes  which  are  associated 
with  enhanced  fluidity. 

Another  aspect  of  the  shear  Induced  structural  distor¬ 
tion  is  presented  by  the  time  dependent  pair  radial  dis¬ 
tribution  function  g(r,  /),  which  is  conveniently  decom¬ 
posed  into  the  following  component  configurational  av¬ 
erages; 

n\r,  t)  =  (r)  +  OvV rl)u„ {r,  t)+  (yl/*f)uy,(r,  t) 

+  (zi/rt)u„(r,  t)  +  Cvz/V))/„(r,  /)  +  •••  .  (11) 

If  4»(r,  t)  is  the  average  of  in  the  radial 

element  r-r  +  dr  about  a  molecule  i  at  time  t,  then 

t)~  15/„(r,  t)/(4n(»*dr)  , 

-  3  Kaa  (r,  t)  +  M»a(r,  t)  +  «^(r,  /)  ,  (12) 

where  ami 

,g« »(r,  t)~  15/aS(r,  t)/ (4 itprzdr)  , 

=  «.*(*-,/),  (13) 

where  a*li. 
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FIG .  9.  The>r(r)  (•'••),  ga  (r>  < — ) .  and  A’,y  (>)  ( - 1  of  I  lie 

sheared  minus  the  unsheared  portions  of  the  n  1.01:101  and 
lxt  =  0.7108  segmented  calculation. 


The /ae(r,  /)  have  been  calculated  at  various  times 
during  the  segment.  The  difference  between  the  gaa  (r,  /) 
of  the  steady  state  sheared  (/)=  1.01304  and  kx,=  0.  7168) 
system  and  the  unsheared  system  at  the  same  density 
bgaa(r,  l)  is  shown  in  Fig.  9.  It  reveals  that  shearing 
produces  a  net  movement  of  particles  in  each  coordina¬ 
tion  shell  towards  the  origin  molecule.  This  is  also  il¬ 
lustrated  by  the  difference  between  the  pair  distribution 
functions  of  the  sheared  and  unsheared  liquids  shown  in 
Fig.  10. 

The  shear  rate  dependent  gx,(r,  /)  at  steady  state  of 
the  p=  0.8442  system  are  compared  with  the  spherically 
averaged  g(r,  t )  from  these  samples  in  Fig.  10.  They 
give  convincing  evidence  of  angular  distortion  in  each 
coordination  shell  such  that  in  the  positive  .vz  quadrants 
the  inside  of  each  shell  is  on  average  depleted  of  par¬ 
ticles  when  compared  with  the  unsheared  liquid.  In  con¬ 
trast,  the  outside  of  each  shell  has  an  excess  of  mole¬ 
cules  in  the  positive  xz  quadrants,  when  compared  with 
the  unsheared  medium.  The  opposite  changes  take  place 
in  the  negative  quadrants. 

The  time  dependence  of  ga,,(r,  /)  and  g„t(r,  t),  shown 
in  Fig.  11,  reveal  that  structural  reorganization  having 
the  symmetry  of  shear  (xz/r2)  is  much  faster  than  that 
of  ( ar/r 2)  symmetry.  Within  the  latter  series,  (x1/^) 
and  (zVx2)  evolve  more  rapidly  to  their  steady  state 
values  than  (y2/r2).  For  example,  when  p=  1. 01304  and 
«„  =  0.  7168,  the  former  three  averages  have  nearly 
reached  their  steady  state  values  by  /-  0. 25  whereas 
that  of  the  (y Vrl>  summation  hardly  differs  from  zero 
by  that  time. 

These  observations  suggest  a  possible  mechanism 
for  the  origin  of  the  Stress  overshoot,  When  a  liquid  is 
sheared,  the  shearing  forces  alter  its  structure  so  that 
the  steady  state  stress  attained  corresponds  to  a  new 
lower  viscosity,  not  that  characterizing  the  starting 
liquid.  At  high  rates  of  shear  the  structural  evolution 
necessary  to  attain  this  new  viscosity  is  significantly 
slower  than  the  viscoelastic  shear  relaxation.  Viewed 
simply  this  means  that  the  original  viscosity  can  re¬ 
main  for  a  sufficiently  long  time  to  enable  the  shear 


stress  to  exceed  its  steady  state  value.  The  situation 
is  more  complicated  than  this  sii  iplo  explanation  pre¬ 
sents,  however,  because  Uie  overshoot  takes  place  at 
times  when  a  purely  viscous  model  for  shear  stress  is 
not  an  entirely  acceptable  approximation. 

Another  probe  into  the  time  dependent  structural 
changes  associated  with  shear  stress  viscoelasticity  is 
provided  by  the  angular  rotation  function  0(r)  defined 
below: 

I  *  *  . 

»(»-)=(EE(VixrIJ)/rj/)AV  ,  (14) 

\  i>j  / 

where  tne  cross  product  is  confined  to  the  .vz  plane. 

When  a  material  is  sheared,  molecules  observe  their 
neighbors  rotate  around  them  with  an  angular  velocity 
such  that  on  average  9{r)~  £„/3.  The  time  and  radial 
dependence  of  tl  ■  function  gives  information  about 
hindered  rotation  and  hence  the  nature  of  structural  dis¬ 
ruption  associated  with  shear  stress  relaxation.  Large 
departures  in  i>(r)  from  ita  continuum  value  of  «„/3  in¬ 
dicate  major  structural  rearrangement  at  that  intermo- 
lecul3r  separation.  Figure  12  shows  b(r)  averaged  over 
three  time  zones  after  the  start  of  s'eady  shear  in  the 
p.-  1.01304,  0.7168  calculation.  This  and  similar 

plots  for  the  calculation  involving  a  step  in  strain  only 
exhibit  a  spike  in  0(r)  at  r  =  0.  S  at  times  when  the  strain 
rate  has  been  turned  off  and  which  is  clearly  above  the 
noise.  This  shows  that  the  slowly  decaying  part  of  the 
stress,  after  a  change  in  applied  strain  rate,  results 
from  further  rearrangements  in  the  inside  of  the  first 
coordination  shell,  which  presumably  require  reorgan¬ 
ization  of  the  liquid  structure  before  they  can  occur. 

The  following  discussion  uses  the  MD  results  to  evalu¬ 
ate  models  for  nonlinear  shear  behavior. 

Previous  MD  studies  of  this  subject  have  used  the 
Hee-Eyrtng  hyperbolic  sine  model2  to  fit  the  shear  rate 
dependence  of  rj.  From  simple  activation  energy  con¬ 
siderations  the  following  relationship  can  be  derived: 

Vt?o=8inh',(c„T)/(c„T)  ,  (15) 

where  rfo  is  the  shear  viscosity  in  the  limit  of  zero 
shear  rate.  The  relaxation  time  r  characterizes  the 


r 


FIG.  10.  The  shear  rate  dependence  ofg(r)  and  g^(r)  for  the 
0.  8-142  states.  t-„.0  (— )i  0.1686  (••••);  0. 3373 

( . )i  e„  =  0.6746  ( - ). 
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HO.  XI.  The  functions  (r,t)  andtf«j<r,  f)  tor  the  calculation  of  Fig.  10:  (a)**,  (b)  yy.  (o)  at,  (d)  x»  for  the  t  =  0. 25  (••••)  and 
(=1.8  ( — )  states. 


medium’s  unsheared  state  and  Is  roughly  an  order  of 
magnitude  larger  than  the  average  shear  relaxation  time 
[=  ft  C(f)df-i7o/G,J.  Although  the  model  relates  r  to  a 
flow  volume  V,  through  the  relationship 

r  =  ijo Vf/2kgT  ,  (16) 

r  is  usually  taken  as  an  adjustable  parameter.  Reason¬ 
able  fits  to  the  MD  results  can  be  achieved  by  taking  Vf 
=  1.1  independent  of  density.  For  example,  for  the  p 
=  1.  01304  state,  this  provides  better  than  15%  agree¬ 
ment  with  the  data;  at  the  other  densities  the  agreement 
ranges  from  13%  to  18%.  The  general  trends  in  the 
plots  of  rj/tjo  versus  «„  at  all  densities  suggest  that  the 
disagreement  is  not  totally  statistical;  the  shear  thinning 
effect  is  systematically  more  pronounced  (i.  e. , 
falls  more  rapidly)  than  Eq.  (15)  would  predict.  Wc  also 
note  that  the  parameter  r  increases  more  rapidly  with 
density  than  does  the  shear  relaxation  time  r(.  At  the 
lowest  density  rf  r, »  20  whereas  for  the  highest  density 
state  t/t,»  100. 

A  perhaps  more  satisfying  model  has  Its  basis  in  a 
trend  in  the  shear  rate  dependence  of  g„(r,  t),  shown  in 
Fig.  10.  The  magnitude  of  the  angular  distortion 
roughly  increases  in  p:  ^portion  to  the  shear  rate,  as 
does  the  shear  stress.  This  suggests  that  some  of  the 


properties  of  a  sheared  liquid  can  be  characterized  by 
a  quasielastic  solid  model  in  which  the  only  nonequilib¬ 
rium  property  is  a  time  dependent  recoverable  shear 


r 

FIQ.  12.  The  normalized  angular  rotation  tunotion  3#(r)/i„ 
averaged  over  the  time  regions  0.17  to  0.33  (— ),  1.64  to  2.06 

< - ),  and  3.86  to  4,48  (••••)  for  tbs  segmented  p>  1.01304 

and  rM  ■  0, 7108  calculation. 
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FIU.  13.  A  comparison  between  the  spherically  averaged  pair 
correlation  function ^(t)  (—1,  gt,{r)  (•  •*•),  and  -  Ae^rlifc/dr) 
( - )  for  the  p=  1 .  and  t  a  =  0.0-183  state. 


strain  A *,,(/).  Such  an  approach2  is  obtained  from  a 
Taylor  expansion  of  g(r,  l): 

uir,  /)  -A-oh  )  +  c  a «„(/)  ) 

+  y  +  •  •  •  .  (17) 

Hence, 


A'(r.  l)  =  gt,{r)  +  xz  A «„(/)av>-  +  x'z'  A«;,(/)(^ 7r 

-i'$/rJ)/2  +  **  ±e‘,(t)fii/(2r)  -*•••,  (18) 

where  jfJ  and  denote  dg^/dr  and  rfVo /rfr*,  respective¬ 

ly.  A  comparison  between  Eqs.  (11)  and  (18)  reveals 
that 

»'A«„(/)/rf  =  «„(r,  f)  .  (19) 

For  most  of  the  density/ shear  rate  combinations  stud¬ 
ied,  a  Ac „(/)  can  always  be  chosen  so  that  Eq.  (19)  is 
an  excellent  approximation.  Good  agreement  for  the  p 
■=  1.  2603  and  «„=  0.  048  calculation  is  illustrated  in 
Fig.  13.  The  fair)  was  approximated  with  little  error 
by  the  spherically  symmetrical  g(r)  of  that  calculation. 
The  recoverable  strain  necessary  to  satisfy  Eq.  (19)  was 
found  to  be  the  Hooke’s  law  strain  (=  o,,/Gm)  to  a  good 
approximation.  Similar  calculations  were  performed 
for  the  other  states.  The  g0(r)  were  taken  from  the  low¬ 
est  shear  rate  at  each  density  or  where  possible  from 
the  unsheared  portion  of  the  segment  In  non-steady  state 
calculations.  Thii  latter  method  was  adopted  in  order 
to  explore  a  larger  region  of  phase  space  than  would  be 
possible  with  a  single  simulation  at  equilibrium,  es¬ 
pecially  at  these  high  densities.  The  best  fits  for  Aa„ 
are  again  given  with  few  exceptions  by  i7„,/G.,  as  Table 
HI  reveals. 

Equations  (17)- (19)  also  predict  changes  in  the  con¬ 
figurational  parts  of  the  total  energy  and  normal  pres¬ 
sure  components  that  take  the  form 


6  £/„(/)--- 


and 


,)  ••  A’o/V)/T5 

o 

+  pl.fAc^(/)  /  S  flr< I, '>■)/, /3+  ■■ 

0 


(20) 


ciAc; 


M 


-W-jA/r) 


where  /!„  3/105,  Ayy  .  1/105;  and  BIX  =  Byy-.-.  1/15, 

1/5. 


1’hr  predictions  for  )  aid  using  this 

model  and  ‘ iking  •  ‘rainn  fr.-m  flq.  My),  arc  several  r>' 
ders  of  magnitude  less  than  those  obtained  by  MD.  It  is 
thought  that  this  is  because  the  model  only  takes  into 
account  those  changes  in  "bulk”  properties  (e.  g. ,  Vc 
and  Paa)  which  arise  directly  from  distortion  of  the 
liquid  having  gtl  symmetry.  Clear  evidence  of  other  in¬ 
duced  strains,  having  gaa  symmetry,  has  already  been 
presented.  Although  these  are  derived  from  the  shear¬ 
ing  action,  they  are  not,  at  present,  related  to  it  in  a 
well  defined  way.  It  is  apparent  that  in  the  MD  results 
the  major  source  of  the  6i’c  and  6Paa  originates  from 
this  latter  structural  alteration  of  the  liquid  which  is  not 
incorporated  in  the  above  model.  The  success  of  this 
theory  for  predicting  shear  stress  follows  from  the  rela¬ 
tively  small  change  in  g(r)  that  occurs  on  shearing  the 
fluids— even  in  the  nonlinear  region.  In  addition,  the 
recoverable  shear  strain  (which  can  be  considered  to  be 
an  expansion  argument)  induced  is  small  (*0.07).  It 
suggests  that  models  characterizing  the  nonlinear  shear 
stress  can  be  constructed  using  a  recoverable  shear 
strain  as  the  only  nonequilibrium  parameter. 


The  analysis  in  the  preceding  paragraph  should  be  re¬ 
garded  as  a  tentative  hypothesis  to  be  explored,  rather 
than  as  a  firm  conclusion.  It  is  not  unreasonable  to  ex¬ 
pect  that  the  shear  induced  changes  6f/„  and  6Paa  will 


TABLE  III.  A  comparison  between 
the  fitted  recoverable  strain  of  Eq. 
(19)  and  the  equivalent  elastic 
medium  strain  om/Gm. 


p 

A<m 

tr„/G. 

U.  8442 

0.0422 

0.005 

0.006 

0.8442 

0.1086 

0.027 

0,019 

0.8442 

0.3373 

0.037 

0.  037 

0.8442 

0.0745 

0.004 

0.064 

0.  92802 

0.0435 

0.01 

0.01 

0.  92862 

0.1741 

0. 1)20 

0.024 

0.  92802 

0.3482 

0.042 

0.042 

0.  92802 

0,  0003 

0.009 

0.069 

1.01304 

0.0800 

0.  024 

0.024 

1.01304 

0.  1792 

0. 034 

0.034 

1.01304 

0.3584 

0.044 

0.044 

t.  01304 

U. 7108 

0.  070 

0.070 

1.2003 

0.0183 

0.041 

0.041 

i .  2003 

0. 1030 

0.U39 

0.03B 

1.2003 

0.3801 

0.046 

0.050 

1.2663 

0.7721 

0,051 

0.069 
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l'ltl.  II.  The  nme  dependence  ol  the  a  *  components  ol  the 

stress  ( — )  and  susceptibility  ( - )  tensors  for  the  p=  1.20(13 

system  subjected  to  the  strain  rate  history  ea- 0.7731,  0  W 
'  3.  cu(«*  0,  2.2<i  i.-l. 


depend  (perhaps  rather  strongly)  on  system  size.  If  this 
is  indeed  the  case,  one  might  find  that  the  discrepancies 
between  the  measured  changes  and  those  calculated 
from  Eqs.  (20)  and  (21)  are  artifacts  resulting  from  the 
small  size  of  the  MD  system. 

Experimentally,  the  stress  produced  by  large  ampli¬ 
tude  strains  can  be  measured  using  the  optical  technique 
of  birefringence. A  fundamental  assumption  of  the 
analysis  requires  that  the  stress  and  susceptibility  ten¬ 
sor  S  respond  in  a  fixed  ratio  to  the  shear  strains.  In 
order  to  test  this  assumption  the  components  of  S, 
i.  e. , 


,  *  *  _  ,  * 
iti  M  l 


4ira! 

“  y 


I»!  *li 


have  been  calculated.  In  Eqs.  (22)  and  (23),  a  is  the 
molecular  polarizability,  which  is  assumed  to  be  equal 
to  o3  here.  A  stress -optical  coefficient  C(  is  defined  to 
be 

C0=:ii„/2o„«,  ,  (24) 

where  n»  is  the  refractive  index  of  the  medium.  The 
stress -optical  coefficient  is  found  to  be  a  l/G.  for  many 
materials.  U,M  in  order  to  test  this  assumption,  a„(f) 
and  Slt(/),  taken  from  the  p=  1. 2083  and  i„=0.  7221 
segmented  calculation,  are  shown  in  Fig.  14.  Even  un¬ 
der  these  highly  nonlinear  conditions  this  proportionality 
is  obeyed  quite  welt  This  is  perhaps  not  surprising 
when  the  similarity  of  the  configurational  average  for 
each  is  considered.  Table  I  reveals  that  Ct  is  relatively 
insensitive  to  shear  rate  when  compared  with  tj,  in 
agreement  with  experimental  evidence,  "  and  to  a  lesser 
extent  density. 

The  susceptibility  is  a  probe  which  weights  longer 
range  Interactions  to  a  greater  extent  than  the  stress. 
This  is  evident  in  comparable  amm(t)  and  Smm(t)  profiles. 


The  ia9(f)  do  not  exhibit  the  same  definite  time  depen¬ 
dence  as  does  the  &’„„(/).  The  latter  is  dependent  on  the 
interaction  potential  and  is  thus  more  sensitive  to  the 
short  range  structural  distortions  accompanying  shear. 
Nevertheless,  these  calculations  give  support  for  the 
use  of  this  technique  in  studying  stresses  in  the  nonlin¬ 
ear  regime. 

For  each  normal  tensorial  property  Tum  (e.g. ,  stress, 
recoverable  strain,  an-,1  susceptibility)  the  direction  and 
its  perpendicular  along  which  the  minimum  and  maxi¬ 
mum  value  of  that  property  occur  define  the  minor  and 
major  axes  of  the  tensortal  distortion  ellipsoid.  The 
smallest  angle  x  between  either  of  these  two  directions 
and  the  v  direction  is  to  a  first  approximation  the  sr.ne 
for  stress,  recoverable  strain,  and  susceptibility.  The 
expressions 

X  =  tan-'(27„/(r„-rM)l  (25) 


X  -  tan*1  (2/ A«„)/2  <26> 

are  readily  derived  and  predict  that  tor  the  states  stud¬ 
ied  \  differs  little  from  45°. 

til.  CONCLUSIONS 

We  have  performed  nonequilibrium  molecular  dynam¬ 
ics  calculations  on  U  fluids  with  the  aim  of  gaining  in¬ 
sights  into  the  mechanism  of  liquid  failure  at  high  levels 
of  shear  stress.  The  shearing  action  has  been  found  to 
change  the  liquid  structure  so  that  there  is  a  tendency  to 
stratify  along  the  lines  of  shear,  enabling  flow  to  take 
place  more  easily.  The  resulting  increase  In  fluidity  is 
evident  in  the  enhanced  self -diffusion  coefficients  and 
the  rate  of  decay  of  the  relaxation  functions  relative  to 
the  unsheared  parent  system.  In  other  words,  the  fluid 
has  a  shorter  memory  of  its  past  behavior. 

The  structural  distortions  giving  rise  to  the  shear 
stress  chaises  are  relatively  small  and  are  to  a  large 
extent  decoupled  from  other  structural  changes.  They 
can  be  characterized  using  a  first  order  perturbation 
expansion  from  the  equilibrium  structure.  The  only 
nonequilibrium  parameter  needed  is  a  recoverable  shear 
■train  which  is  given  to  a  good  approximation  by  the 
shear  stress  divided  by  the  shear  rigidity  modulus.  Un¬ 
fortunately,  we  do  not  understand  how  to  formulate  the 
mechanism  by  which  the  shearing  changes  the  other 
aspects  of  the  liquid’s  structure,  which  undoubtedly  has 
a  strong  influence  in  determining  the  recoverable  shear 
strain.  Model  predictions  of  “bulk"  property  changes, 
e.g. ,  in  configurational  and  normal  pressure  compo¬ 
nents,  are  similarly  only  at  an  elementary  stage  of  de¬ 
velopment. 
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Viscoelastic  properties  of  a  L-  -inard-.lones 
i*  (LJ)  model  system  studied  in  a  previous  work 

were  compared  with  the  behavior  of  inorganic 
and  organic  qiass-forminq  liquids  and  solid 
glasses.  In  the  region  of  small  stresses  and 
strains  both  the  LJ  model  system  and  "real 
materials"  demonstrate  linear  behavior. 
After  adjustment  of  time  scales  the  stress 
relaxation  curves  for  Ld  model  and  fused 
>  silica  were  found  to  be  similar  despite  tne 

16  orders  of  magnitude  difference  in 
viscosities  of  the  liquids. 

At  hiqh  stresses  the  viscosity  of  the  L.l 
model  system  drops,  stress  relaxation  occurs 
faster,  and  a  peak  appears  in  the  stress 
versus  time  plots  for  shearing  with  a 
constant  rate.  This  phenomenon,  well-known 
for  hiqh  polymers,  has  also  been  observed  for 
inorqanic  glasses.  The  curves,  &  log  n 
vs.  log  o,  that  show  a  drop  of  viscosity  of 
Rb  O  -  Sio  glass  and  LJ  model  at  stresses 
higher  than  10*  Pa,  are  in  remarkable 
agreement.  The  LJ  model  is  shown  to  provide 
a  semi-quantitative  reprcsentat ion  of  a 
variety  of  glass-forming  liquids. 


INTRODUCTION 

Inorganic  glass-formi nq  liquids  are  often  subjected  to  severe 
manufacturing  and  operating  conditions.  For  example,  during  forming 
operations  glass  melts  experience  shear  strains  and  rates  about  10* 
cm/cm  and  10*  sec-1  ,  respectively  —  enormous  values.  Calculations 
of  melt  behavior,  however,  are  normally  made  using  experimental  data 
obtained  at  strains  and  rates  10  orders  of  magnitude  lower.  On 
lookinq  for  methodR  to  study  the  behavior  glass-forming  liquids 
subjected  to  high  amplitude  disturbances  wo  came  to  a  molecular 
dynamics  (MD)  investigation  of  a  model  liquid  system. 
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MO  i3  a  computer  simulation  technique  for  prol)in<)  the  mlcrtvuv'-j;  • 
behavior  of  a  system  that  allows  one  to  account  for  the  many-ho.tv 
nature  of  the  molecu  ar  interactions.  it  involves  to! lowing  so: 
representative  number  of  molecules,  n,  by  solving  the  classical 
equation*  of  motion  on  a  high  speed  computer  to  obtain  their 
trajectories  in  time  steps  at  -  10“*".  From  the  results,  of  an  Mb 

experiment  —  essentially  a  record  of  the  system's  path  thtough 
phase  space  --  one  can  determine  any  measurable  property  of  the 
system,  as  well  as  many  that  are  experimentally  inaccessible. 

Heyes,  Kim,  Montrose  and  l.ltovit*  |l|  studied  by  this  method  the 
model  system  of  108  Leonard- Jones  (id)  particles.  Those  "molecules" 
interact  through  pairwise  additive  forces  described  hy  the  IJ  fi-1,: 
potential 

♦  ( r )  -  -4u  I (a/r)k  -  (a/r) 1  *  I ,  (1) 

where  r  is  the  separation  of  an  interacting  pair  and  u  and  t  are 
constants  with  dimensions  of  energy  and  distance  respectively.  The 
method  of  computer  simulation  was  shown  IP  to  be  a  very  effective 
technique  for  atudytnq  the  behavior  of  an  IJ  model.  In  this  paper 
we  compare  the  results  obtained  for  the  LJ  system  with  the  behavior 
of  inorganic  and  organic  qlaas-forming  liquids  to  assess  the  extent 
to  which  it  can  be  considered  representative  of  real  materials. 
Although  no  real  qlass-formlnq  liquid  can  he  characterized  hy  a 
potential  of  the  LJ  form  (only  the  noble  gas  fluids  Ar,  Kr,  etc,, 
are  quantitatively  represented  in  this  way)  there  may  be  some 
important  common  features  in  dynamical  behavior  of  liquids  which  are 
essentially  independent  on  the  chemical  compos i t ion  and  the  extent 
of  supercool inq  relative  to  the  molting  point. 

• 

SHEAR  STRESS  RELAXATION,  LINEAR  BEHAVIOR 

Consider  first  the  linear  {small  stress  and  strain)  behavior. 
Figure  1  shows  the  stress  response  when  a  shear  strain  is  switched 
on  at  the  Instant  t  *  0  and  is  kept  constant. 


Figure  i 

A  Comparison  of  the  Normalized  Stress  Relaxation 
Functions  for  the  LJ  Model  System  and  Sl.>2  mass 
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Figure  2 

The  Time  Dependence  of 
the  Reduced  Stress, 

0  ( t )/£  ,  for  the  LJ  model: 
_  I  -  4.4  .  10* ®  . 

- f  -  1.66  «  10“  s“>  , 

....  I  ■  3.32  «  10' 1  •-> . 

VISCOSITY I  NON— LINEAR  BEHAVIOR 


STRESS  RELAXATION  NON-LINEAR 
BEHAVIOR 


Experiments  on  stress  relaxation 
for  constant  strain  r..  Imposod  at 
t  »  0  wore  carrion  out  as  a 

function  of  1  .  There  is  little 
change  in  the  "instantaneous  shear 
rigidity  aa  a  function  of  t  > 
howcvcr(  there  is  a  marked 
difference  In  the  relaxation 
function,  t(t).  The  relaxation 
function  in  the  hiqh  initial 
strain  cate  decays  more  rapidly 
to  toro  at  long  times.  Precisely 
tho  tame  concluaions  were  drawn 
by  the  authora  (61  who  atudied 
extensively  the  stress  relaxation 
in  poly  isobutylene.  Their  data 
demonstrate  the  dependence  of  the 
atoraqe  modulus  C  on  radian 
frequency  w  for  different  maximum 
amplitudes  of  deformation  rate 

*max  ’  fo*  .  v.  .ls  the 

frequency).  Increase  of  (or 

c0  for  any  qiven  v)  cuts  of*  the 
lonq  time  tall  of  G'(w>  function 
similarly  to  that  observed  in 
this  work  for  ♦ (t) . 


Piqure  3  shows  experimental  data  for  4  systems  with  viscosities 
ranqlnq  from  10”1*  to  10' *•  Pa*  a.  Despite  an  enormous  difference  in 
structures  for  all  systems  we  have  similar  behavior)  that  is,  a 
decrease  of  viscosity  beqinninq  at  a  certain  ranqe  of  values  of  the 
stress.  The  stresses  at  which  deviations  beqin  are  different  and 
clearly  dependent  on  the  details  of  structure.  In  this  connection 
it  is  interesting  to  note  the  coincidence  of  the  curves  for 
Inorganic  glass  and  the  LJ  model. 

STRESSES  IN  A  LIQUID  SHEARED  WITH  A  CONSTANT  RATE  i 


Typical  curves  obtained  for  the  I.J  model  at  small  and  large  t  ure 
shown  In  Pig.  2.  The  fact  that  at  high  i  the  stress  versus  time 
curve  rises  to  a  maximum  value  and  then  decays  to  its  lonq  time 
limltlnq  value  is  well  known  for  real  materials.  The  curves,  quite 
similar  to  ones  shown  In  Fiq.  2,  were  reported  for  poly(vinyl 
acetate)  |9|,  and  other  qlass-forminq  liquids. 


DISCUSSION 

1.  Experimental  evidence  is  provided  to  support  the  assumption  that 
there  are  important  basic  features  in  viscoelastic  behavior  which 
are  common  to  materials  with  enormous  differences  ir.  structure. 
These  are:  a)  transition  from  elastic  behavior  at  t  ►  0  to  viscous 
flow  at  larqer  times.  The  transition  itself  can  be  characterized  by 
fairly  similar  distributions  ot  relaxation  times  for  LJ  model 
compared  to  that  of  real  materials:  inorganic  qlasses,  inorqanic  and 
organic  monomeric  liquids.  h)  in  all  systems  studied  similar 
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Adjusting  time  scales  we  compare  relaxation  of  stresses  (in  portions 
of  initial  atresa)  in  fused  silica  glass  ami  the  Ut  model.  For  SIO 
we  used  the  data  of  Leko  and  Mescheryakova  |2|.  Taking  into  account 
the  16  orders  of  magnitude  difference  in  viscosities  the  similarity 
in  behavior  seems  to  be  remarkable.  The  UT  curve  indicates  that  the 
stress  relaxation  in  the  !«?  model  in  characterised  by  the 
distribution  of  relaxation  times.  The  U3  curve  in  Pin.  1  is  well 
approximated  by  a  function  of  the  double  exponential  form 

G(t)  ■  G.|  (1  -  a)  expl-t/tj  |  ♦  a  expl-t/t  f  1|  ,  (2) 

where  G*  is  the  instantaneous  shear  rigidity  and  a  the  constant. 
The  ratio  tj/t.,  which  in  usually  considered  an  a  characteristic  of 
the  distribution  width,  is  equsl  to  10.  This  ratio  ir  very  close  to 
that  found  in  most  inorqanie  and  organic  monomeric  liquids  (1,4). 
For  example,  for  B  03  ,  •  9,  and  for  boroaillcate  glass,  i  /r 


Coincidence  of  the  parameters  describing  macroscopic  behavior  does 
not  mean  the  identity  of  microscopic  mechanisms.  It  means,  however, 
that  the  varieties  of  mechanisms  arc  similar  in  the  sense  that  much 
faatar  mechanisms  coexist  with  slower  ones. 

VISCOSITY i  LINEAR  BEIIAV’OR 

The  next  step  is  to  apply  to  the  U)  model  some  of  the  equations  of 
the  linear  thaory  of  viscoelasticity  which  were  proven  to  be  valid 
for  Inorganic  glasses.  The  equation  relatinq  the  normalised  shear 
stress  relaxation  function,  4(t),  the  shear  rigidity  modulus,  ,  to 
the  viscosity  of  the  material,  is 

» 

n  “  Q-  J"  dt  ♦  C  t ) .  (1) 


In  Raf.  (1|  a  was  computed  for  the  UT  system  in  this  fashion,  l.a., 
using  GU  and  *(t)  obtained  in  the  stress  relaxation  experiment  at 
constant  strain  c .  Then  the  experiment  on  UJ  model  was  carried  out 
to  determine  the  stress  response  of  the  system  when  a  constant  shear 
rate  I  is  initiated  at  t  ■  0. 

The  stress  was  found  to  bshave  as  shown  in  Fig.  2  (see  the  upper 
curve).  Again  using  the  linear  theory,  this  response  is  described 
by  the  form 

o(t)  -  l  Q.  J*  dt«  4 (t *  > ,  (4) 

which  gives  Eq.  (3)  at  t  ■  -i  8(t*»)/<  •  n  •  G.  /"  dt  #(t).  The 
initial  slope  is  G»i  and  the  lonq-time  value  1$  n.  c.  The  shear 
rigidity  and  viscosity  determined  in  this  second  type  of  experiment 
agreed  to  within  atatistical  uncertainty  with  that  found  in  the 
first  type  of  experiment.  This  procedure  or  similar  ones  were 
carried  out  for  many  commercial  silicate  qlasses  1 5 J  to  show  that 
they  behave  aa  linear  viscoelastic  bodies  at  low  stress  and  strain. 
Thus  the  UT  model  is  quite  representative  of  many  real  materials, 
the  behavior  of  which  at  low  awlttudes,  is  described  by  linear 
theory. 
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non-linear  effects  were  observed  at  higher  amplitudes  of  shearinq 

strain  or  stress. 


logo  (Pa) 


Figure  3 

Viscosity  as  a  Function  of  Stress  for  Several  Liquids  and  Glass,  a 

(the  data  for  Rb  O-SiO  glass  was  obtained  by  Li  and  Uhlmann  [7] 
the  data  for  polymer  solutions  are  qiven  by  Ferry  in  (8|). 

2.  The  drop  of  liquid  viscosity  with  increasing  stress  or  strain 
rate  seems  to  be  a  general  phenomenon;  it  occurs  in  all  liquids 
studied  under  shear,  uniaxial  extension  and  compression,  and  in  some 
liquids  even  under  hydrostatic  pressure. 

All  viscoelastic  functions  and  parameters  as  defined  refer  to 
Isothermal  changes  of  state.  At  the  Bame  time  tne  enerqy, 

dissipated  during  viscous  flow  may  heat  the  sample.  The  work  of 

viscous  shearing  heating  is  W  ■  oc  «  oft  »  o2t/n,  which  indicates 
that  the  dissipated  energy  increases  as  a2.  Therefore,  the  control 
of  constant  temperature  of  the  sample  is  an  Important  feature  of  the 
measurements  of  viscosity  as  a  function  of  stress.  While  it  is  a 

difficult  problem  in  some  "real”  experiments  is  was  fairly  easy  to 

keep  constant  the  temperature  of  the  "sample"  during  MD  experiments. 
Thus,  the  factor  which  might  be  considered  as  a  reason  for  the  drop 
in  viscosity  has  been  excluded.  L-rge  deformations  are  likely  to 
produce  distortion  or  reorganization  of  the  system's  structure. 
Evidence  supporting  this  hypothesis  was  obtained  in  Ref.  |lj  by 
defining  and  measuring  a  set  of  planar  distribution  functions  for 
the  LJ  model  system.  The  representation  of  the  modified  structure 
i 8  sketched  in  Fig.  4  (neglect  at  first  the  line  connecting  the 
center  of  the  molecules). 
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Prom  this  one  can  see  the  nature 
of  the  structural  reorganization 
that  has  occurred t  the  system  has 
arranged  itself  into  layers 
aligned  with  t  lie  flow  direction. 
The  connection  of  the  center  of 
the  molecules  by  the  line  gives  an 
impression  of  the  analogy  in  the 
behavior  of  simple  liquids  and 
polymers  of  high  molecular  weight. 
Figure  4b  shows  the  influence  of 
shearing  upon  the  two  polymer 
chains  with  entanglement  coupling. 
The  nature  of  stress  or  strain 
induced  structural  rearrangements 
is  still  quite  speculative. 
Figure  4  gives  one  of  the  possible 
pictorial  interpretations  of  the 
similarity  in  behavior  which  has 
been  found  in  this  work  for  tin- 
materials  with  different  s tinc¬ 
tures. 
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A  Pictorial  Representation 
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Deformations  in  Simple  Liq¬ 
uids  (a)  and  Polymers  with 
Entanglement  Coupling  ;h). 
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ABSTRACT 

A  molecular  dynamics  investigation  of  the  microscopic 
behavior  of  simple  liquids  subjected  to  shearing  conditions 
similar  to  those  found  in  elastohydrodynamic  rolling  contacts  is 
outlined.  The  behavior  of  a  model  liquid  and  more  complicated 
traction  fluids  at  high  levels  of  stress  are  strikingly  similar. 
The  calculations  give  insights  into  the  structural  origin  of 
shear  failure  and  reveal  a  restructuring  of  the  compacted 
material  into  layers,  which  significantly  .  reduces  its  shear 
viscosity.  This  change  is  also  manifest  in  an  increase  of  bulk 
pressure  and  decrease  in  high  frequency  shear  rigidity  modulus 
and  thermal  conductivity.  Thermal  effects  in  the  nonlinear 
region  of  viscosity  are  consistent  with  standard  macroscopic 
theory . 


can  be  done  by  appending  the  idea  of  a  limiting  shear  stress  to 
the  low  amplitude,  i , '  linear,  shear  response.  In  order  to 
test  this  hypothesis  and,  more  generally,  to  provide  a  framework 
on  which  to  base  analytical  representations  of  lubricant 
viscoelastic  behavior,  we  have  undertaken  a  Molecular  Dynamics 
(MD)  investigation  of  a  model  liquid  system  under  con  :  cions 
approximating  those  found  in  an  END  contact.  Even  though  the 
model  system  is  particularly  simple,  one  can  expect  to  reproduce 
the  major  features  of  shear  characteristics  of  real  lubricants. 
Then,  because  an  MD  experiment  provides  one  with  a  complete 
microscopic  record  of  the  system's  evolution  in  time,  one  can 
"measure11  aspects  of  its  behavior  that  are  not  accessible  in 
conventional  experiments,  but  which  provide  useful  insights  into 
the  origin  of  the  properties  of  liquids  under  the  above  extreme 
conditions. 

The  MD  calculation  procedure  follows  rather  closely  that 
pursued  in  certain  model  calculations  of  the  contact  [3]  in  that 
we  observe  a  portion  of  liquid  as  it  is  sheared.  This  is 
equivalent  to  concentrating  on  a  "slice”  of  lubricant  as  it 
traverses  the  region  of  contact.  Although  a  simultaneously 
applied  time  dependent  pressure  and  shear  rate  could  have  been 
studied,  we  have  confined  our  attention  at  this  stage  to 
considering  the  response  of  a  model  fluid  to  a  time  varying 
shear  rate  at  a  series  of  densities  (pressures). 


INTRODUCTION 


In  several  recent  papers  Bair  and  Winer  (1,2)  have  examined 
the  problem  of  understanding  the  shear  stress  response  of 
liquids  in  concentrated  contacts  with  particular  attention  to 
the  possibility  of  predicting  elastohydrodynamic  (BHD)  traction 
data.  They  point  out  that,  under  the  severe  conditions  to  which 
the  lubricant  is  subjected  — ■■  large  shearing  rates  and  large, 
rapidly  varing  normal  pressures,  there  is  a  dearth  of  primary 
laboratory  data  on  which  to  base  a  physical  model  of  the 
lubricant's  behavior.  In  a  series  of  well  conceived 
experiments,  they  have  attempted  to  remedy,  at  least  partially, 
this  situation.  Using  their  data,  they  were  then  able  to 
develop  a  physically  appealing  phenomenological  model  of  the 
non-linear  shear  response  of  liquids  that  is  reasonably 
successful  in  predicting  EHD  traction  results. 

In  order  to  predict  traction,  one  must  be  able  to  compute 
the  average  shear  stress  in  the  contact,  which  is  equivalent  to 
knowing  the  time  evolution  of  the  stress  in  a  fluid  element  as 
it  moves  through  the  contact  zone,  where  it  encounters  a 
simultaneously  applied  steady  shear  rate  and  a  time  varying 
pressure.  The  response  of  the  system  under  these  conditions 
can,  for  low  shearing  rates,  be  parameterized  in  terms  of  the 
shear  viscosity,  n*  a  shear  rigidity  modulus,  Q. ,  and  a 
distribution  of  relaxation  times.  At  the  rather  large  shearing 
rates  generally  encountered,  this  relatively  simple  description 
must  be  modified  to  incorporate  the  non-linear  character  of  the 
dynamical  stress  response.  Bair  and  Winer  (2]  suppose  that  this 


TUB  MODEL 


MD  is  a  computer  simulation  technique  for  probing  the 
microscopic  behavior  of  a  system  that  allows  one  to  account  for 
the  many-body  nature  of  the  molecular  interactions.  It  involves 
following  some  representative  number  of  molecules,  N,  by  solving 
the  classical  equations  of  motion  on  a  high  speed  computer  to 
obtain  their  trajectories  in  time  steps  of  typically  0.01  ps 
[4] .  Since  only  a  rather  a  limited  number  of  molecules  is 
treated  (generally  less  than  a  thousand),  the  molecules  are 
confined  to  remain  in  a  cubic  box  which  is  surrounded  by  images 
of  itself  to  avoid  severe  boundary  effects.  In  many 
circumstances  the  artificialities  introduced  by  this  procedure 
ate  rather  insignificant.  From  the  results  of  an  MD 
"experiment"  —  essentially  a  record  of  the  system's  path 
through  phase  space  —  one  can  determine  any  measurable  property 
of  the  system,  as  well  as  many  that  are  experimentally 
inaccessible. 

The  model  system  under  investigation  is  an  assembly  of  108 
Lennard-Jones  (LJ)  particles;  these  molecules  interact  through 
pairwise  additive  forces  described  by  the  LJ  6-12  potential, 

♦  (r)  -  -4u  [ (a/r)6  -  (a/r)12]  ,  (1) 

where  r  is  the  separation  of  an  interacting  pair,  and  u  and  a 
are  constants  with  dimensions  of  energy  and  distance, 
respectively.  All  the  computed  quantities  are  in  so-called  IJ 
reduced  units  which  are  given  in  terms  of  u,  a,  and  m,  the  mass 


of  an  individual  particle.  Although  no  real  lubricant  can  be 
characterized  by  a  potential  of  the  LJ  form  (only  the  noble  gas 
fluids  Ar,  Kr,  etc.,  are  quantatively  represented  in  this  way) 
their  structural  and  dynamical  behavior  can  be  expected  to  be 
qualitatively  similar  to  that  of  the  model  fluid.  As  a  result 
one  can  expect  that  the  microscopic  insiqhts  derived  from  MD 
experiments  will  serve  as  useful  guidelines  for  understanding 
and  modeling  the  behavior  of  real  lubricants  and  perhaps  for  the 
"molecular  engineering"  of  improved  liquid  traction  fluids  and 
lubricants.  For  the  purpose  of  comparison  with  experiment 
frequent  conversion  from  reduced  to  "real"  units  will  be  made, 
using  the  parameters  for  Ar  which  are  given  in  Table  1. 

Two  methods  of  shearing  the  box  of  molecules  were  adopted. 
The  results  of  each  method  agree  within  statistical  error; 
however,  each  method  is  well  suited  to  examine  a  particular 
aspect  of  the  investigated  phenomenon. 

The  first  method  attempts  to  mimic  the  environment  a 
lubricant  experiences  in  EHD  rolling  contacts,  in  which  the 
fluid  forms  a  thin,  almost  parallel  film  of  length  400^m  and 
thickness  lym,  typically  under  fully  flooded  conditions.  This 


thin  film 

is 

achieved  in 

the  model  system 

by  employing 

periodicity 

in 

the  x  and  y 

directions 

only. 

so  that  an 

infinitely 

thin 

film  is  simulated  IS).  The 

moving  rollers  are 

replaced  by 

two 

fluid  layers 

translat ing 

with 

the  desired 

velocities,  Uj  and  U2 .  Molecules  from  the  three  regions  are  not 
allowed  to  mix  and  are  kept  from  doing  so  by  reflection  boundary 
conditions  in  the  z  direction,  which  are  discussed  in  detail 


-*>- 


elsewhere  [5].  The  fluid  walls  (FW)  are  maintained  at  a  preset 
temperature  by  adjusting  their  velocities  at  each  time  step  of 
the  integration  scheme.  The  main  features  of  the  model  are 
pictorially  represented  in  Fig.  1. 

The  second  method  is  a  modification  of  the  Homogeneous 
Shear  technique  (4]  ( MHS ) ,  which  is  better  suited  to  study  the 
viscoelastic  behavior  of  the  fluid  because  it  employs 
periodicity  in  all  directions  so  that  a  bulk  material  is 
sheared.  It  also  enables  a  desired  strain  rate  to  be 
instantaneously  achieved  throughout  the  liquid.  The  molecular 
trajectories  are  disturbed  from  those  governed  by  equilibrium 
fluctuations  by  imposing  a  shear  strain  rate,  e ,  on  the  system. 
This  is  accomplished  by  displacing  the  x-coordinate  of  the  i'th 
molecule  in  the  MD  cell  by  an  amount  4x^  »  where  z^  is 
the  a  coordinate  of  particle  1  and  at  is  the  length  of  the  time- 
step,  for  as  long  as  the  shearing  is  required.  The  results  of 
this  technique,  in  particular  those  relating  to  viscoelasticity 
(derived  using  a  time  dependent  strain  rate),  complement  those 
of  the  FW  method. 

We  have  used  non-equi librium  MD  experiments  to  investigate 
the  linear  and  non-linear  shear  response  of  the  LJ  system  under 
a  variety  of  thermodynamic  conditions.  The  temperature,  T,  and 
number  density,  p,  initially  studied  were  close  to  the  normal 
freezing  temperature,  Tfl ,  of  the  LJ  liquid,  that  is,  0.722  kBT/u 
and  0.8442  Na3/V,  respectively.  Densif ications  of  10%,  20%,  and 
50%  relative  to  the  starting  density  were  achieved  by  suddenly 
(within  one  time  step)  compressing  the  system  at  constant 
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temperature.  The  densif ications  were  achieved  so  rapidly  that 
the  systems  were  kineticaily  prevented  from  crystallizing.  The 
reduced  number  densities  of  the  other  states  considered  are  thus 
0.92862,  1.01304  and  1.2663. 

Each  calculation  was  undertaken  for  a  duration  of  at  least 
30  reduced  time  units  (65  ps).  The  computations  were  carried 
out  on  a  PDP  11/40  computer. 
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RESULTS  AND  DISCUSSION 

We  first  consider  the  t  dependence  of  3hear  viscosity  on 
shear  rate  in  the  non-linear  region. 

In  the  MHS  method  a  uniform  x  velocity  profile  was  set  up 
instantaneously  whereas  a  period  of  several  time  units  was 
required  for  the  fluid  walls  to  drive  the  central  region  into 
this  state.  An  x  velocity  gradient  was  achieved  by  maintaining 
the  upper  and  lower  walls  at  velocities  of  U  and  -U, 
respectively.  It  is  relevant  here  to  mention  that  at  all 
densities  and  shear  rates  considered,  the  fluid  film  sheared 
uniformly  in  the  z  direction  without  slip  at  the  walls,  and  not 
concentrated  in  the  center  as  has  been  speculated  before  [6] . 
An  example  of  a  vx(z)  so  derived, from  the  p  ■  0.8442  and  c  = 

0.3373  calculation  is  shown  in  Pig.  2.  The  value  of  e 
corresponds  to  an  extremely  large  shear  rate  —  on  the  order  of 
lO^a**1  depending  somewhat  upon  the  particular  values  of  u,  a, 
and  m  that  are  chosen.  Such  larqe  values  of  £  were  chosen  since 
it  is  desirable  to  explore  the  system's  behavior  at  shear  rates 
on  the  order  of  the  reciprocal  of  the  viscoelastic  relaxation 
times.  For  the  LJ  fluids  these  are  approximately  unity  (about 
10“12  s).  The  shear  stress  needed  to  define  a  shear  viscosity 
is  readily  derived  from  the  velocities  and  positions  of  the 
molecules.  The  a0  component  of  the  stress  tensor  is  defined 
below. 
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where  V  is  the  volume  of  the  MD  cell,  is  the  a  velocity 
component  of  molecule  i  relative  to  the  average  velocity  and  a  ^ j 
is  the  a  component  of  the  separation  between  molecules  i  and  j. 
The  subscripts  for  the  shear  stress  component,  o xz ,  are  dropped 
in  further  discussion. 

At  low  shear  rates  (AU/h)  plots  of  traction  versus  slip  are 
linear,  however  at  high  stress  levels  (••  C^/100)  the  traction 
increases  less  rapidly  than  slip  speed.  In  other  words  the 
shear  viscosity  (■  o/c )  decreases  from  the  equilibrium  value  of 
n0  with  inreasing  shear  rate.  The  MHS  method  was  used  to  obtain 
n  versus  c  for  a  variety  of  densities.  The  forms  of  n/nQ  versus 
stress  in  Pig.  3  resemble  closely  those  obtained  from  twin-disk 
experiments  for  the  fluids  5P4E  (7J  ,  L63/1271  [8)  and  Oxilube 
8 .S/140  [9]  after  making  the  appropriate  conversions  to  real 
uni'-s.  The  ratio  n/n0  descends  markedly  from  unity  in  the 
region  **  50MPa.  The  de  cease  in  n  with  c  is  not  due  to  thermal 
heating  because  the  calculations  were  conducted  isothermally . 
Another  method  of  presenting  the  behavior,  logfn/n0)  vs. 
loglex],  where  t  is  a  characteristic  shear  relaxation  time  (» 
n q /G» ) ,  again  shows  marked  quantitative  similarity  between  the 
model  fluid  and  a  real  fluid,  which  is  in  this  case  5P4E  ( 2 ) . 
These  results  suggest  that  the  model  and  real  liquids  might 

i 

share  a  common  mechanism  of  faii  ire  at  high  levels  of  shear 
stress.  With  the  support  of  the  previous  results  we  now  discuss 
associated  changes  in  the  fluids  which  until  now  have  only  been 
speculated  upon. 
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These  high  stresses  cause  extreme  distortions  of  the 
system's  structure.  The  intermolecular  forces  are  not 
sufficiently  strong  to  support  the  enormous  stress  levels  that 
develop}  consequently  there  is  a  rupturing  and  reorganisation  of 
the  structure  resulting  in  enhanced  flow  behavior.  Evidence 
supporting  this  hypothesis  was  obtained  by  measuring  the  number 
density  in  the  s  direction  between  the  fluid  walls.  The  density 
profile  of  Fig.  5,  derived  for  p  ■  1.01304  and  t  ■  0.3584 
suggests  the  form  of  this  structural  reordering.  The  molecules 
align  themselves  into  "glide"  planes  along  the  line  of  shear.  A 
pictorial  representation  of  the  structure  suggested  is  given  in 
Fig.  6.  This  layering  was  not  observed  when  the  walls  were 
stationary  or  moving  in  the  same  direction  with  equal  velocity. 
In  this  context  it  is  relevant  to  note  that  IR  studies  of  flow 
under  EHD  conditions  have  also  been  interpreted  in  terms  of 
molecular  alignatent  along  the  line  of  shear  [10]  . 

These  structural  rearrangements  are  also  manifest  in 
changes  in  the  so-called  normal  pressure  components  ,  Poa  (» 
-o00 ) .  Typically,  Psz  >  PXJt  >  Pyy  although  their  differences 
are  small  when  compared  with  each  component's  change  on 
shearing.  The  ss  pressure  component  increases  are  shown  in 
Fig.  7  and  predict  an  improvement  in  the  load  carrying  capacity 
equivalent  to  ~0.5GPa  at  most. 

It  is  now  widely  accepted  that  knowledge  of  the 
viscoelastic  response  of  lubricants  is  necessary  to  interpret 
EHD  results  over  the  range  of  possible  operating  conditions  [9] . 
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In  order  to  better  understand  non-linear  viscoelasticity  the  MD 
calculations  were  conducted* with  a  time  varying  strain  rate. 

We  consider  the  dynamical  responae  of  the  L3  material  to  an 
imposed  shear  that  is  "switched  on"  suddenly  using  the  MHS 
method.  Two  types  of  MD  experiments  are  of  interest!  (a)  the 
response  of  the  system  to  a  steady  shearing  rate  (uniform 
velocity  gradient)  £  •  i vx/a i  that  is  initiated  at  time  t  ■  0* 

and  (b)  the  response  of  the  system  to  a  constant  shear  strain  (a 
pulse  of  shearing  rate)  imposed  at  time  t  «  t'.  Linear  response 
theory  provides  us  with  relations  between  the  two  types  of 
responses!  the  shear  stress  response  to  a  step  of  shear  strain 
can  be  written  as 

• 

o(t)/c  «  C?^<t)#  (3) 

where  4(t)  is  the  normalized  (4(0)  ■  1)  shear  stress  relaxation 
function.  In  terms  of  these  parameters  the  shear  stress 
response  to  a  "small"  steady  shear  rate  imposed  at  t  *  0  is, 

a(t)/£  -  G.  /J  dt*  4(f).  (4) 

since  the  steady  state  (t  ♦  • )  value  of  this  is  just  the 
viscosity  of  the  material,  n0 ,  we  have, 

n0  «  G „  JJ  dt*  4(f)  -  G.t  .  (5) 


The  last  equality  serves  to  define  the  shear  relaxation  time  t. 


Experiments  of  the  type  described  as  (a)  above  were  then 
carried  out  as  a  function, of  e  for  values  of  e  well  beyond  the 
range  for  which  linear  behavior  can  be  expected.  A 
representative  selection  of  the  results  is  displayed  in  Figs.  8 
and  9  in  which  the  density  (constant  c  )  and  shear  rate  (constant 
density)  dependence  of  n  are  shown,  respectively. 

Several  features  of  these  figures  are  worthy  of  comment. 
At  short  times  the  system  behaves  largely  as  an  elastic  solid  in 
that  the  stress  increases  linearly  with  time, 

o(t)  -  G„et.  (6) 

The  shear  rigidity  modulus  rises  linearly  with  pressure  at 

♦ 

roughly  the  same  rate  as  more  complicated  fluids,  as  is  revealed 
in  Fig.  10.  Liquids  cannot  permanently  support  a  shear  strain 
and  consequently  the  stress  levels  off  to  a  value  determined  by 
the  shear  viscosity  of  the  sheared  state.  There  is  evident 
"shear  thinning,"  or  a  decrease  of  viscosity  with  increasing 
shear  rate.  The  predicted  value  of  the  stress  build-up  for  the 
p  *  0.92862,  e  “  0.6963  state  using  equation  (4),  shown  in 
Fig.  8,  clearly  overestimates  the  actual  stress  obtained. 

Also,  a  maximum  in  the  stress  versus  time  curve  for  those 
simulations  in  the  high  density/shear  rate  regime  is  observed. 
This  indicates  that  the  large  shear  deformations  are  producing 
structural  distortions  which  lag  in  time  somewhat  behind  the 
growth  of  the  shear  stress  but  follow  rather  closely  the 
observed  time  variation  of  the  normal  pressure  components  in  the 
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material  (see  Fig.  11).  This  relatively  slower  development  of 
what  can  be  thought  of  as  transient  stress-induced  glide  planes 
is  consistent  with  the  observed  shear  thinning  and  stress 
"overshoot"  behavior  presented.  It  is  also  consistent  with  the 
observation  that  the  stress  relaxation  function  decays  more 
rapidly  in  a  highly  sheared  system  than  in  an  equilibrium 
system.  A  comparison  of  the  two  types  of  behavior  is  presented 
in  Fig.  12.  The  form  of  these  4>(t)  consists  of  a  rapidly  (t  < 
0.2)  decaying  portion  which  is  presumably  inertial  in  origin  and 
due  tc  small  motions  of  the  molecules  within  the  random  network 
structure.  Further  stress  relaxation  requires  the  cooperative 
rearrangement  of  larger  molecular  groups  and  is  consequently 
slower.  The  attenuation  of  * (t)  with  increasing  shear  rate  is 
indicative  of  the  enhanced  stress  relieving  properties  of  the 
structurally  reorganized  system.  The  stress  overshoot  observed 
in  non-linear  viscoelastic  experiments  on  polymers  [11]  has  a 
form  that  is  quite  similar  to  the  curves  obtained  for  the  LJ 
system.  While  the  microscopic  mechanisms  for  the  two  systems 
are  obviously  quite  different,  the  general  principle  that 
associates  the  stress  overshoot  phenomena  (as  well  as  shear 
thinning)  with  a  time  dependent  structural  reorganization 
process  remains  quite  valid. 

In  this  context  it  is  perhaps  relevant  to  note  that  a  more 
rapid  decay  of  the  stress  relaxation  function  can  be  thought  of 
as  both  a  suppression  of  the  long  relaxation  times  associated 
with  stress  relaxation  and  also  a  reduction  in  the  shear 
rigidity  modulus,  at  low  frequencies,  G,  (a  rather  well  known 
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rthenomenon  in  polymeric  systems  (12)).  It  is  found 

experimentally  that  the  extracted  shear  rigidity  modulus  of  the 

* 

liquids  confined  in  concentrated  contacts  is  up  to  an  order  of 
magnitude  smaller  than  those  observed  obtained  from  Independent 
laboratory  (ultrasonic)  measurements  at  the  same  temperature  and 
pressure.  The  large  shear  rates  (~  1 0**  s-1  )  that  exist  in  the 
inlet  zone,  even  under  "pure  rolling"  conditions,  could  be 
sufficient  to  reorder  the  liquid  preferentially  along  the  line 
of  shear.  The  observed  increase  in  shear  modulus  with  rolling 
speed  [13]  for  a  wide  range  of  fluids  is  particularly 
significant  evidence  for  this  interpretation  because  if 
molecular  orientation  of  the  molecule  along  the  line  of  flow  is 
the  cause  of  this  phenomenon  then  it  would  be  expected  that  the 
modulus  would  increase  with  rolling  speed  because  there  is  less 
time  for  molecular  rearrangement  to  take  place. 

Thermal  effects  are  a  major  factor  in  limiting  the  ability 
of  a  fluid  to  sustain  a  shear  stress  at  high  strain  rates  [6] . 
The  PW  technique  is  well  suited  to  determine  the  validity  of 
using  standard  macroscopic  theory  on  thin  films.  Prom 
conservation  of  energy  considerations  the  temperature  profile 
under  conditions  of  steady  shear  flow  as  given  by  [14] , 

K  *l2±l±  -  -  n  e 2 /(p  C ) ,  (7) 

3  Z2 

where  K  *  A/(pC),  and  A  and  C  are  the  thermal  conductivity  and 
specific  heat,  respectively.  Assuminq  a  temperature  independent 
shear  viscosity  then  the  solution  to  equation  (7)  is 
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T ( z )  -  T(0 )  -  n  e2  z2/(,2X  ).  (8) 

A  typical  temperature  profile,  taken  from  the  p  *  1,01304 
and  c  -  0.3584  state,  is  shown  in  Pig.  13.  All  the  temperature 
profiles  can  be  fitted  rather  well  to  the  above  form  which  is 
parabolic  in  z  and  hence  can  be  used  to  estimate  X .  Note  that 
these  macroscopic  equations  only  apply  to  a  steady  state  viscous 
material  and  do  not  apply  under  viscoelastic  conditions  where 
some  energy  is  stored  as  potential  energy  and  is  not  dissipated. 
The  density  and  shear  rate  dependence  of  X  is  given  in  Pig.  14. 
The  trends  are  consistent  with  the  general  behavior  that  thermal 
conductivity  decreases  as  a  material  becomes  more  fluid-like. 


CONCLUSIONS 


We  have  performed  non-equilibrium  molecular  dynamics 

t 

calculations  on  LJ  fluids  with  the  aim  of  gaining  insights  into 
the  mechanisms  of  liquid  failure  at  high  levels  of  stress  such 
as  is  found  in  BHD  contacts. 

The  model  liquids  "rupture"  at  similar  stresses  to  those  of 
more  complicated  real  fluids  due  to  marked  structural  reordering 
into  layers.  Paradoxically,  although  structural  aspects  of  the 
liquids  present  a  more  solid-like  appearance,  other  more 
dynamically  related  properties  such  as  shear  viscosity  and 
thermal  conductivity  manifest  changes  which  are  associated  with 
enhanced  fluidity. 
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Table  1 

A  summary  of  che  reduced  units  in  terms  of  which  all  quantities  are 
given.  Boltzmann's  constant  is  denoted  by  kR. 


Quantity 

Reduced  Unit 

ij 

31  Unit  for  Ar 

distance 

a 

0.3405  nm 

mass 

m 

6.64  x  10“26  kg 

energy 

u 

1.65  x  10~21  J 

time 

a(m/u ) 1 ' 2 

2.16  ps 

density 

a"3 

42.1  kg  mol/m3 

temperature 

u/kB 

119.8  K 

pressure,  stress 

u/a3 

41.8  MPa 

modulus 

u/a3 

41.8  MPa 

viscosity 

(mu)1  /2/a2. 

9.d3  x  10~5  Pa  s 

thermal  conductivity  kB(m/u)"1 /2o~2 

1.88  x  IQ-*2  J  nT1  K~ 
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Captions  to  the  Figures 

Fig.  1  A  pictorial  representation  of  the  relationship  between 
the  FW  MD  technique  and  an  EHD  contact.  The  encounter 
of  of  molecules  with  the  reflection  boundarier  in  the  z 
direction  and  periodic  boundaries  in  the  x  and  y 
directions  are  illustrated. 

Fig.  2  The  z  dependence  of  the  average  x  velocity  for  the  p  = 
0.8442  and  c  =  0.3373  state  using  the  FW  model.  The  z 
coordinate  system  is  measured  from  the  center  of  the 
film. 

Fig.  3  The  non-linear  behavior  of:  A  BP  L63/1271  [8]  (P  *  0.8 
GPa,  T  -  27'C,  U  -  1.12  m/s),  0  5P4E  f 7 1  (P  *  0.45  GPa, 
T  -  27*C,  U  -  0.6  m/s)  and  9  Oxilube  85/140  [9]  (P  «  1.2 
GPa,  T  •  30°C,  U  ■  2.2  m/s)  from  twin-disk  experiments. 
The  MHS  models  A  p  *  0;8442,  7  p  *  0.92862,  Op* 

1.01304  and  open  square  p  -  1.2663. 

Fig.  4  The  non-linear  behavior  of:  5P4E  (2)  r  *  9  x  10~3  s  from 
experiment.  The  filled  in;  i  p  *  0.8442,  V  p  «  0.92862, 
0  p  «  1.01304  and  square  p  »  1.2663  from  MHS  MD.  The 
relaxation  times,  t,  are:  0.12,  0.30,  0.43,  and  1.62, 
respectively. 

Fig.  5  The  density  profile  from  the  p  -  1.01304  and  e  -  0.3584 

FW  calculation. 

Fig.  6  A  pictorial  representation  of  the  structural  changes 
that  take  place  on  going  from  an  unsheared  (a)  to  a 
sheared  (b)  state.  The  arrows  denote  the  line  and  plane 
of  shear. 


Fig.  7  The  variation  in  the  zz  normal  pressure  conponent  for 
the  MHS  systems.  A  p  =  0.8442,  V  p  ■  0.92862,  0  p  • 

1.01304  and  square  p  »  1.2663. 

Fig.  8  The  time  dependence  of  the  reduced  stress,  a  (t)/c  i  — 
p  *  0.8442,  e0  *  0.6745;  »•»»  (a)  actual  ••••  (b) 

predicted  using  equation  (4)  p  ■  0.92862,  c0  ■  0.6963; 

- P  -  1.01304,  c  Q  *  0.7168;  -  p  »  1.2663,  - 

0.7721.  The  associated  strain  rate  variations  are 
indicated  above  the  figure. 

Fig.  9  The  time  dependence  of  the  reduced  stress,  o(t)/c0,  for 

the  p  »  1.01304  state:  -  eQ  ■  0.0896,  — —  tQ  » 

0.3584,  *•••  cfl  =  0.7168. 

Fig.  10  The  pressure  dependence*  of  :  A  Santotrac  40  19] ;  O 
di ( 2-ethylhexyl )  phthalate  |9]  from  a  twin-disk 

apparatus.  Solid  line  »  (0.293  +  1.70P)  GPa 

di( 2-e) phthalate  [15],  The  filled  in  O  are  obtained 


from  MHS  MD. 


Fig.  11  The  time  dependence  of  the  normal  pressure  components 
for  the  p  ■  0.92862  and  eQ  »  0.6963  state:  ..... 

»  ^xj]  ( t ) ,  and  — — — —  Pyy(t). 

Fig.  12  The  shear  stress  relaxation  functions  for  the  p  ■  0.8442 

c  ■  0  and  -  c  ■  0.6745.  The  step  in  strain  (c  ■ 

3.6  *  10“3 )  at  t  ■  0  is  illustrated  above  this  figure. 
Fig.  13  The  variation  in  temperature  in  the  z  direction  for  the 
p  ■  1.01304  and  e  ~  0.3584  state  where  TQ  ■  0.722. 
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14  The  effect  of  non-linearity  on  thermal  conductivity  (X). 
The  symbols  refer  to  the  same  calculations  aa  for 
Fiq.  7. 
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FIGURE  13 


Section  IV.  A  theory  of  non~linear  response  In  liquids  end  Amorphous 
solids. 


by  C.  J.  Montrose 
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A  THEORY  OF  NONLINEAR  RESPONSE  IN  LIQUIDS 
AND  AMORPHOUS  SOLIDS  * 


by 


C.J.  Montrose  ♦  > 

Vitreous  State  laboratory  -#4 
The  Catholic  University  of  Ameri.u 
Washington,  D.C.  20064 

ABSTRACT 

A  mathematical  model  of  nonlinear  viscoelastic  stress  relaxation  is  derived 
using  heuristic  arguments  relating  shear  response  to  dynamical  strain  Induced 
structural  changes  in  the  material.  A  parameter  characterising  the  limiting 
steady-state  shear  stress  appears  as  a  central  feature  of  the  model.  General 
formulae  are  presented  and  some  representative  specific  cases  are  considered.  The 
predicted  behavior  is  in  qualitative  agreement  with  both  MD  and  "real" 
experimental  observations. 


*  This  work  was  supported  by  the  Office  of  Naval  Research  (Contract  No. 
N00014- 81-K-02 96 ) . 


The  shear  stress  response  at  time  t  to  a  ehear  strain  rate  history  imposed 
over  the  time  Interval  0  to  t  is  given  for  small  stresses  hy  the  linear  convolution 
integral  form 

o(t)  -  |dt'  G0(t-t’)fc<t').  (1) 

In  this  o(t)  is  the  shear  stress  at  time  t,  e(t')  ■  dc/dt'  is  the  shear  strain 

rate,  and  G0(t-t')  is  the  shear  stress  relaxation  function.  In  terms  of 

this  last,  one  can  obtain  the  usual  shear  viscoelastic  parameters  of  the  material: 

(1)  The  Instantaneous  shear  rigidity  modulus  is  just  the  value  of  Go  at 
the  zero-value  of  its  argument,  i.e. 

G0(O)  -  G„  .  (2) 

We  make  use  of  this  to  define  a  normalized  shear  stress  relaxation 
function  gfl(s)  -  - 

Gq (s)  5  0,80(8)  .  (3) 

(2)  The  shear  viscosity  is  the  integral  of  the  relaxation  function  G0 : 

oo 

n0  ■  ^ds  G0(s)  (4) 

(3)  The  shear  relaxation  time  is  the  integral  of  the  normalized  relaxation 
function: 

T0  -  Jds  g0(s)  -  no/G«  »  (5) 

the  second  equality  following  from  the  use  of  Eq.  (3)  in  the  first. 

That  the  relations  expressed  in  Eqs.  (2)  -  (5)  are  appropriate  can  be  seen  by 
considering  a  feu  special  strain-rate  histories. 

Consider  first  a  step  of  strain  AeQ  switched  on  at  the  instant  of  time  -  0+. 

The  shear  strain  rate  is  e(t’)  “  Aeofi(t’)  and  from  (1)  the  stress  respons  is 

o(t)  -  AcqGq (t) .  (6) 

At  the  instant  t*0,  this  gives  the  instantaneous  elastic  response.  Thus,  Go(0) 
is  indeed  the  Instantaneous  shear  rigidity,  G,. 


2 


In  a  similar  fashion,  we  may  examine  the  response  to  a  steady  strain  rate 
£q  switched  on  at  time  =  0+;  here  the  stress  is  given  by 

t  t 

a(t)  -  e0^dt’  Gq (t-t * )  «  e0£ds  Go(s),  (7) 

the  second  equality  following  from  substituting  t'**t-s  in  the  first.  The 
steady-state  viscous  behavior  is  found  by  taking  the  t  -*■  “  limit: 

00 

a(t  -*  °°)  *  eQ|ds  G0(s)  v'8) 

The  ratio  o(t  «®)/£0  is  the  usual  definition  of  the  viscosity  in  agreement 
with  Eq.  (4). 

The  results  and  formulae  ju^t  given  are  the  usual  viscoelastic  equations 
for  the  linear  response  of  a  system  when  the  stresses  and  strains  (and  strain-rates) 
are  small.  The  purpose  of  this  paper  is  to  develop  formulae  analogous  to  these 
which  can  be  used  in  situations  of  large  stresses  and  strains,  and  which  will 
of  course  reduce  to  these  in  the  appropriate  limits.  The  paper  consists  of  two 
main  parts:  in  the  first  a  heuristic  approach  is  used  to  derive  a  generalized 
stress  response  formula  analagous  to  (1)  and  to  use  this  to  define  effective,  i.e. 
strain-rate  dependent,  moduli,  viscosities,  and  relaxation  times.  The  behavior 
for  some  particularly  simple  impressed  strain  rate  histories  is  examined.  In 
the  second  part,  specific  forms  (based  roughly  on  the  MD  results  obtained  in  this 
laboratory1)  of  the  relevant  shear  and  structural  response  functions  in  the 
theory  are  assumed.  The  consequences  of  these  are  worked  out  and  the  resulting 
behavior  is  compared  with  some  existing  experimental  work. 

mm 

In  generalizing  Eq.  (1)  for  the  case  of  finite  amplitude  shear  rates,  we 
assume  that  the  same  form  can  be  used  if  explicit  allowance  for  the  dependence  of 
(and  dynamic  response  of)  the  stress  relaxation  function  on  strain  rate  is  made. 


Specifically,  It  is  assumed  that  the  stress  relaxation  function  depends  upon  the 
structural  configuration  of  the  material,  which,  in  turn,  responds  to  the  impressed 
shearing  rate.  The  shear  stress  is  thus  expressed  as 

t 

cr(t)  =  /dt0G[t-t  |t,e]e(t0).  (9) 

0 

The  form  of  G  for  a  giver,  shear  rate  is  taken  to  be 

G [ t~t0  j  t , e]  =  Gq ( t — t o ) {  1  -  X/dtjf  (tj)F^t-ti  |t,e] }.  (10) 

0 

where  F  is  the  structural  relaxation  function  (explicitly  time  shear  rate 
dependent)  and  is  a  "strength  of  coupling"  parameter,  for  which  a  more  physical 
interpretation  will  be  given  below.  The  structural  relaxation  function  is 
assumed  to  depend  on  the  material  configuration  in  the  same  manner  as  G,  that  is, 
we  take 

F  [t— t !  1 1 ,  tj  =  fo(t-t1)U  -  X/dt2e(t2)F[t-t2!t,eJ  }.  (11) 

0 

Combining  (9)- (11)  gives  for  the  shear  stress 

t  t  t 

o(t)  =  /dt0e(t0)G  (t-t0)  -  x/dt0e(t0)G0(t-t0)  /dt1e(t1)  f Q(t-t1) 

0  0  O' 

t  t  t 

+  X2/dt0e(t0)G0(t-t0)/dt,e(t1)f0(t-t1)/dt2e(t2)F(t-t2| t,s)  (12) 

0  0  0 

Repeated  use  of  Eq .  (11)  in  (12)  leads  to  an  infinite  series  form  ior  the  stress: 
t  °°  t 

o(t)  =  /dt0e(t0)G0(t-t0)  I(-X)n|/dt,c(t’)f0(t-t')|n  (13) 

0  n=0  0  J 

It  should  be  appreciated  that  the  form  in  (13)  can  be  obtained  more  formally 

from  the  general  Green-Rivlin  integral  equation,2 

t  t  t 

o(t)  =  /dt0e(t0)C0(t-t0)  +  /dt0e(t0)/dt14(t,)G.  (t-t  ,  t-t.) 

0  0  o  ' 

t  t  t 

+  /dt0t(t0)/dt1e(t1)/dt2e(t2)G2(t-t0,  t-tj,  t-t2)  +  .  (14) 

0  0  0 
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by  assuming  that  the  n^1  order  response  function  is  fnc.torizable  as 
Gn(t-t0,  ...  t~tn)  =  G0{t-t0)(-X)nf0(t-t1)f0(t-t2)  ...£0(t-tn). 


(15) 


With  the  usual  restrictions  on  the  summand,  the  series  in  (13)  can  be 
summed  in  closed  form*,  the  result  is 


o(t)  - 


/dt0E(t0)G0(t-t  ) 
0 


(16) 


i  +  x/dt'e(t')i0(t-t’) 

0 


Eq .  (16)  is  the  fundamental  dynamical  result  of  this  paper.  Note  that  without 
any  loss  of  generality  we  may  assume  that  -^o  Is  normalized  as 


j0(O)  -  1, 


(17) 


absorbing  any  other  factors  into  the  coupling  parameter  X.  This  also  implies 
a  definition  of  the  structural  relaxation  time  ts  as 


(18) 


Eq.  (16)  can  now  be  written  in  terms  of  somewhat  more  physical  quantities  by 
examining  the  steady-state  (t  -*■  «)  form  predicted  by  this  equation  when  a  steady 
shear  rate 


0  (t<0) 

e  (t>0) 


(19) 


is  applied.  The  result  is 


<">  '  e0t0Gw/(1  +  c.Xtg) 


(20) 


For  large  shear  rates  this  becomes 


o(») 


t  a- 


^  Gwt q/X'Cq  —  o  l 


(21) 


o*  is  seen  to  be  the  maximum  steady  state  shear  stress  that  can  be  sustained  by 
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the  system.  Eq.  (16)  can  now  be  cast  as 

t 

G„/dt0E(t0)g0(t-t0) 

d(t)  *  - ® - £ - —  (22) 

1  +  t0(G./o*)I  /dt'Kt'jfoft-t') 

Tso  * 

By  considering  the  strain  rate  history  given  in  (19)  it  is  possible  to  define  an 
effective  shear-rate-dependent  viscosity 

n«0)  s  o(~)/e0  ;  (23) 

the  form  of  this  is  easily  found  by  substituting  the  second  equality  in  (21) 
into  (20): 

_  (24) 

n<eo>  *  1  +  £t0  (Gja*) 

This  is  sketched  in  Figure  1.  An  effective  rigidity  modulus  can  also  be 
defined  by  considering  the  instantaneous  stress  developed  upon  application  of 
a  step  in  strain  Ae0»  Here  we  obtain 

G(fe,>  *  -212H - & - — ■  .  (25) 

Acq  1  +  Ae0(G«T0/<r  rs) 


In  the  sample  computations  presented  in  this  section  we  have  used  the 
following  typical  functional  forma  and  parameter  values: 

A.  go (t)  *  a  exp(-t/Ti)  +  (l-a)exp (-t/T2) 
a  -  0.9 
t2/ti  -  9 
t0  -  t2/5 
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B.  j'pU)  =■'  t'xp(-t/r2) 

%  =  T2 

C.  Coo/o*  *  10 

Figure  2  presents  computations  using  these  functions  and  parameters  in  Eq.  (22) 
for  the  strain  rate  history  given  in  (19).  Observe  that  both  the  stress  overshoot 
and  shear  thinning  effects  observed  by  Hoyes  jet  al1  (section  I  of  this  report) 
are  reproduced.  As  expected,  these  effects  become  more  pronounced  for  large 
values  of  the  reduced  strain  rate. 

We  have  also  examined  the  behavior  when  a  step  of  strain  AeQ  is  switched 
on  at  time  =  0.  For  this  rase  the  strain  dependent  shear  rigidity  is  given  by 
Eq .  (25)  and  is  plotted  in  Figure  3.  The  time  evolution  of  the  stress  is 
easily  computed  to  be 


_ _ 8p(f-)  __ 

o(t)  =  G„Ac0  l  +  AE0(GeDT0/o*Ta)^0(t) 


This  is  plotted  in  Figure  4. 

As  a  last  case  we  consider  the  application  of  a  steady  shear  rate  [see  Eq.(19)J 
at  time  =  0,  followed  by  the  application  of  a  small  (infinitesimal)  step  of  strain 
6r  imposed  at  a  time  t*>>  tq,  i.e.  after  the  stress  has  reached  its  steady  state 
value.  The  incremental  stress  So  is  calculated  to  be 


6o( 


G„«e  j  —  r 

C)  =  1  +‘  e0T0(n«/o*)  j80(t)  -  !  +  f qTq (Geo/ o*)  (To/Ts>io<fc)J;  (27> 


in  this  equation,  the  time  variable  t  represents  the  time  following  the 
application  of  the  strain  increment.  From  this  equation  we  can  define  an  effective 
shear-rate  dependent  rigidity 


6o( 0)  1  4  e0T0(G„/o*)(l-Tp/T8) 

<5e  |_1  +  e0^(G„/o*)j2 


G(£0)  H 


(28) 
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and  a  shear-rate  dependent  stress  relaxation  function 

g0(t)  -  eoT0(C,/g*)[g0(t)  -  (T0/T8)i0(t)1 
g(t;feo)  ra  1  +  e0t0(G./(T*)d-T0/T8)  ~  U9' 

This  last  is  plotted  in  Figure  5.  Again  there  is  qualitative  agreement  with  the 
observations  of  Heyes  jRt  al.  * 


1.  D.M. 

2.  See, 
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FIGURE  1 


0.001  0.01  0.1 
NORMALIZED  STRAIN  RATE. 


The  shear-rate  dependent  viscosity  plotted 
versus  reduced  strain  rate.  The  data  points 
are  MD  results  (Section  I,  this  report)  and 
the  value  G/o(limit),  i.e.  G  /a*,  of  13  was 
chosen  to  best  fit  these  data. 
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A  MOLECULAR  DYNAMIC  STUDY  OF  STRUCTURAL  CHANGES 

* 

ACCOMPANYING  MATERIAL  FAILURE t  PRELIMINARY  RESULTS 

by 

S.  D.  Grant  and  C.  J.  Montrose 

Vitreous  State  Laboratory 
The  Catholic  University  of  America 
Washington,  D.C.  20064 

INTRODUCTION 

In  a  recent  paper1  we  reported  a  molecular  dynamics  (MD) 

study  of  the  nonlinear  shear  response  of  stable  and  metastable 

(undercooled)  liquid  systems.  Among  the  principal  results  of 

this  work  was  the  observation  that  accompanying,  and  indeed 

% 

responsible  for,  the  reduction  in  the  effective  viscosity  (shear 
thinning)  at  elevated  shear  rates  is  a  reorganization  of  the 
liquid  structure  into  laminae  approximately  aligned  with  the 
shear  force.  The  effect  of  this  is  to  facilitate  flow  and  hence 
reduce  the  shearing  stress  that  is  required  to  maintain  the  shear 
rate.  At  the  time  it  was  hypothesized  that  this  structural 
reorganization  could  be  viewed  as  the  precursor  of  a  material 
failure  process,  the  actual  fracture  of  the  system  being  prevented 
by  rapid  molecular  diffusion  that  acts  as  a  local  "healing" 
mechanism.  To  investigate  this  hypothesis  we  have  undertaken 
a  set  of  MD  experiments  on  systems  that  have  been  rapidly  quenched 
to  low  temperatures  —  glasses--so  that,  on  the  time  scale  of  the 
experiments,  molecular  diffusion  is  arrested.  In  this  note,  we 
report  some  preliminary  results  of  those  experiments. 


-2- 

I 

I  COMPUTATIONS  AND  RESULTS 

The  system  under  investigation  is  an  assembly  of  108  parti¬ 
cles  (a  few  runs  have  been  carried  out  on  larger  systems)  inter- 

i  j 

acting  pairwise  via  a  Lennard-Jones  6-12  potential 
<t>  (r)  =  4e[  (o/r)  12  -  (o/r)M  . 

!  1 

The  equations  of  motion  were  integrated  using  the  Verlet  algorithm 

2 

with  a  time  step  of  about  0.01.  The  results  that  we  present  here 
,  were  obtained  at  a  temperature  of  0.10  and  a  density  of  1.013. 

^  This  state  was  reached  by  rapidly  (in  one  time  step)  "crushing"  j 

a  triple-point  liquid  system  (T  =  .722  and  p  *  .8442)  to  achieve 
a  ?0%  densif ication,  then  allowing  this  state  to  equilibrate,  and 
►  finally  "quenching"  the  system  (again  in  one  time  step)  to  the 

desired  temperature.  The  rapid  densif ication  and  cooling  steps  1 

I  : 

j  inhibit  those  structural  rearrangements  that  are  required  for 

I  crystallization  to  occur.  1 

|  The  experiment  consists  on  switching  on  a  shear  strain  rate 

b.  at  time  -  03  and  examining  the  subsequent  time  evolution  of  1 

xy  » 

the  system.  The  temperature  was  held  constant  by  scaling  the 
particle  velocities  after  each  time  step.  Among  the  system 

V 

properties  "measured"  were  the  shear  and  normal  stress  components  ■ 

and  distribution  functions  characterizing  the  local  structure,  ? 

Typical  of  the  results  obtained  is  the  plot  of  shea  (xy) 
stress  versus  time  shown  in  Fig.  1.  As  is  evident  the  stress 
rises  until  approximately  t  *  12  at  which  time  (when  the  strain 

is  about  6%)  it  drops  precipitously  to  zero.  This  stress  release  j 

i 

t 

S 
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is  interpreted  as  a  failure  mode  of  the  material;  fracture  in 

the  normal  sense  cannot  occur  because  of  the  confinement 

condition  imposed  by  the  density  constraint  on  the  system. 

This  effect  was  quite  reproducible  occuring  in  three  separate 

runs  with  e  =  .005  and  in  one  with  e  =  .001.  In  all  cases 
xy  xy 

the  stress  release  occurs  at  stresses  «  0.6  and  at  strains  of 
5%  to  6%. 

One  of  the  runs  with  C  -  .005  was  continued  out  to  t  «  75 

xy 

(Fie.  2).  Observe  that  following  the  failure  at  t  s  12,  the 
stress,  after  hovering  near  zero  for  a  short  time  begins  to  rise 
at  t  a  20  before  dropping  at  t  =  42.  Note  that  at  this  point 
the  stress  has  risen  to  a  value  more  than  twice  as  large  as  the 
level  initially  required  to  initiate  failure.  This  growth  and 
sudden  fall  of  the  stress  is  repeated  again  between  t  *  50  and 
t  =  72. 

The  time  dependence  of  tne  other  stress  components  are 
shown  for  the  same  run  in  Figs.  3  through  8. 

Accompanying  each  of  the  major  stress  changes  are  significant 
structural  rearrangements.  These  are  pictured  in  Figs  9(a)  and  (b) . 
It  is  evident  that  the  initial  stress  buildup  causes  a  strati¬ 
fication  of  the  amorphous  structure;  the  failure  of  the  system 
coincides  with  a  reorientation  of  these  strata.  Continued 
shearing  of  the  system  leads  to  a  "reconstruction"  of  the  material 
in  which  it  is  reconfigured  as  a  crystalline  solid. 

This  is  illustrated  in  Figs.  10  -  12  where  the  radial  distri¬ 
bution  function  is  given  at  various  time  instants.  At  the  later 
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time  instants  the  crystalline  character  of  the  material  is  clear; 
the  positions  of  the  first  four  peaks  in  the  RDF's  are  in  the 
ratios  1,  /?,  /J ,  and  /T  indicating  an  FCC  crystal  structure. 
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The  xx~component  of  the  pressure  tensor  (the  negative  of 
the  stress  tensor)  versus  time  following  the  application 
of  a  steady  xy  shear  rate  at  t  -  0. 
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ABSTRACT 

The  viscosity  of  a  soda-lime  silica  glass  was  measured  at  \ 

high  strain  rates.  The  data  show  non-Newtonian  viscous  flow  in 
this  inorganic  oxide  glass  witli  the  viscosity  values  below  the 
expected  Newtonian  value.  Following  the  imposition  of  large, 
steady  strain  rates,  the  observed  stress  increases  with  time  to 
a  maximum  and  then  decreases  to  a  time-independent  value.  A 
comparison  of  the  viscosity  behavior  of  this  glass  with  the 
molecular  dynamics  results  in  a  "Lennard-Jones"  glass  shows  a 
number  of  points  of  correspondence  and  suggests  the  interpreta¬ 
tion  of  the  non-Newtonian  behavior  as  resulting  from  structural 
rearrangements  in  the  material.  The  combined  data  show  that  the 
sustained,  steady-state  stress  asymptotically  approaches  a  maxi¬ 
mum  at  very  high  strain  rates.  This  limiting  stress  is  inter¬ 
preted  as  the  actual  cohesive  strength  of  the  material  and  is 
calculated  to  be  1.4  x  108N/m2  (20,000  psi)  for  the  glass  under 
study . 
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INTRODUCTION 

The  strain  rate  response  of  non-crystalline  materials  under 
an  applied  stress  has  generally  been  observed  to  be  linear 
(Newtonian)  for  low  stress  or  strain  rates.  Organic  materials 
have  exhibited  large  deviations  from  Newtonian  behavior  'Oth  in 
the  pseudoplastic  direction  (below  Newtonian  viscosity)  and  the 
dilatant  direction  (above  Newtonian  viscosity).^  Since  a  cor¬ 
respondence  between  specific  structural  changes  and  the  non- 
Newtonian  response  has  not  been  determined  for  these  materials, 
it  has  been  generally  surmised  that  the  non-Newtonian  behavior 
is  a  result  of  complex  molecular  chain  kinetics  such  as  unfold¬ 
ing,  stretching,  cross-linking,  etc. 

Recent  "experiments"  with  Lennard-Jones  spheres  in  a  glassy 
state  using  molecular  dynamics  (MD)  calculations  have  also  shown 
the  onset  of  non-Newtonian  behavior  under  increasing  strain 
rates  or  applied  stresses. The  calculations  show  that  when 
an  applied  shear  strain  rate  is  imposed  on  a  box  containing 
Lennard-Jones  spheres,  the  resulting  stress  is  proportional  to 
the  applied  strain  rate,  c,  only  for  low  strain  rate  values 
(2t0<<1,  where  t„  is  the  shear  isothermal  relaxation  time  in  the 
Newtonian  region).  At  higher  strain  rates,  the  measured 
stresses  and  the  corresponding  viscosity  do  not  reach  their 
expected  Newtonian  value,  with  the  deviation  growing  for 
increasing  strain  rates.  The  resulting  behavior  is  shown  in 
Figure  1.  Since  in  Lennard-Jones  glasses  the  atoms  interact 
only  via  central  forces,  it  seems  clear  that  the  observed  non¬ 
linear  behavior  in  MD  calculations  is  a  fundamental  property  of 
the  liquid  or  glassy  state. 
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Oxide  glasses  consist:  of  ions  which  bind  with  both  an  ionic 
and  a  covalent  character.  Silicate  glasses  have  not  demon¬ 
strated  any  chain  characteristics,  and  therefore,  they  offer  a 
set  of  materials  whose  behavior  can  be  compared  to  ND  Lennard- 
Jones  glasses.  Some  differences  are  expected  since  silicate 
glasses  exhibit  some  covalent  (directional)  bonding.  But  much 
similarity  is  also  expected  since  a  large  portion  of  their 
molecular  make-up  is  ionic. 

Measurements  by  Li  and  Uhlmann  on  rubidium  silicate 
glasses*  have  shown  the  existence  of  a  non-Newtonian  region  in 
the  viscoelastic  response  of  the  glass  to  an  applied  constant 
load.  As  is  the  case  for  MD  glasses,  the  deviation  from  Newto¬ 
nian  behavior  was  in  the  pseudoplastic  direction.  However,  the 
dynamics  of  the  non-Newtonian  »behavior  of  the  glass  were  not 
studied  in  detail.  We  have  analyzed  their  data,  and  a  compari¬ 
son  of  their  results  to  MD  experiments  and  to  our  results  is 
presented  below. 

In  this  paper,  a  stable  silicate  glass  was  selected  for  an 
investigation  of  non-Newtonian  behavior.  The  glass  is  a  stan¬ 
dard  reference  material  for  viscosity  issued  by  the  National 
Bureau  of  Standards.^  Its  Newtonian  viscosity  has  been  measured 
at  great  length  and  its  stability  in  air  at  high-temperatures 
has  been  established  through  numerous  tests. 

The  experiment  which  we  present  is  conducted  differently 
from  that  of  Li  and  Uhlmann,  and  is  designed  to  reproduce  the 
conditions  of  the  MD  calculations.  The  experiment  is  conducted 
at  constant  strain  rate  rather  than  constant  stress  or  applied 
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load.  While  these  conditions  are  interchangeable  when  the 
material  is  exhibiting  Newtonian  behavior,  they  are  not  equiva¬ 
lent  when  the  response  is  non-linear.  Constant  strain-rate 
experiments  were  performed  because  they  allow  a  detailed  study 
of  the  non-linear  behavior  before  the  onset  of  material  failure. 
Experiment 

Measurements  were  conducted  on  an  oxide  glass  fiber  whose 
composition  is  shown  in  Table  1.  Fibers  were  made  by  hand  draw¬ 
ing  from  the  molten  glass  at  1150°C,  and  had  diameters  of  0.1  to 
1mm. 


TABLE  1*  -  GLASS  COMPOSITION 


Si02 

70.5% 

sb2o3 

1.1% 

K.O 

7.7% 

•  S03 

0.2% 

Na20 

8.7% 

^2^3  1  Fe2^3 

0.2% 

CaO 

11.6% 

*This  glass  is  known  as  NBS-710  Viscosity  Standard. 


Each  tested  fiber  was  inserted  through  a  furnace  12cm  long  with 
a  narrow  central  channel  and  was  attached  to  a  stationary  fix¬ 
ture  at  the  bottom  end  and  to  a  moveable  load  cell  at  the  top 
end.  See  Fig.  2  for  a  schematic  of  the  fiber  and  furnace 
arrangement.  The  furnace  was  designed  to  yield  a  relatively 
constant  temperature  over  most  of  its  length  (±3°C  over  90mm) 
and  to  drop  rapidly  in  temperature  at  each  end  (AT  *  -20°C  in 
less  than  6mm).  This  guaranteed  a  constant,  well  defined  hot 
zone.  At  selected  temperatures,  the  fiber  was  elongated  at  a 
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constant  rate  while  the  resulting  force  was  measured.  The 
measurements  were  conducted  at  all  the  different  strain  rates 
studied  before  the  temperature  was  changed. 

Since  the  length  of  hot  fiber  was  fixed  by  the  furnace, 
this  resulted  in  a  constant  applied  strain  rate.  The  fiber 
diameter  was  reduced  as  its  total  length  increased,  therefore  it 
was  necessary  to  calculate  a  correction  for  the  steady-state 
cross-sectional  area  reduction.  Steady-state  refers  to  the  con¬ 
dition  where  the  cross-sectional  area  decreases  uniformly  over 
the  entire  hot-length  of  the  fiber.  This  is  selected  to  differ¬ 
entiate  it  from  localized  necking  which  we  sha-ii  discuss  later. 
An  exact  solution  to  the  steady-state  correctiw:.  is: 


A(t)  «  AQe"e  *t-to  * , 


(1) 


where  AQ  is  the  area  at  tQ  and  c  is  the  strain  rate.  The  vali¬ 
dity  of  this  correction  was  verified  by  the  measurement  of  the 
expected  Newtonian  viscosity  at  low  strain  rates  before  and 
after  substantial  elongation  at  high  strain  rates.  Using  this 
area  correction,  plots  of  the  variation  of  developed  stress  with 
time  were  obtained  for  different  applied  constant  strain  rates. 
A  typical  dependence  of  stress  on  strain  rate  is  shown  in  Fig. 
3.  Note  the  linear  region  at  low  strain  rates,  and  the  pseudo¬ 
plastic  deviation  at  higher  strain  rates.  The  apparent 
viscosity  n  was  calculated  as  follows; 
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where  a  is  the  developed  stress.*  Figure  4  shows  the  temperature 
dependence  of  the  Newtonian  viscosity  measured  at  low  strain 
rates. 


♦This  is  a  commonly  used  definition  of  viscosity  in  the  Newto¬ 
nian  region.  The  definition  of  viscosity  in  non-Newtonian 
regions  becomes  somewhat  ambiguous  and  others*  have  suggested 
that  n«l/3(Jo/8c)  is  a  more  appropriate  definition  (the  two 
definitions  are  equal  in  the  Newtonian  region).  The  use  of  Eq. 
2,  however  allows  us  to  plot  viscosity  as  a  function  of  time  at 
fixed  i  and  this  is  more  useful  for  our  fixed  strain  rate  exper¬ 
iments.  The  use  of  the  differential  form  requires  data  with 
very  low  scatter  to  calculate  the  viscosity  with  any  accuracy 
and  therefore  is  not  used  here. 
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Discussion 

A.  Experimental  Results 

A  study  of  the  non-linssr  behavior  of  this  glass  and  a  com¬ 
parison  with  the  Lennard-Jones,  molecular  dynamics  experiments 
are  best  effected  by  looking  at  the  time  dependence  of  the 
apparent  viscosity  for  a  given  applied  constant  strain  rate. 
Figures  5  and  6  show  the  non-Newtonian  behavior  for  two  differ¬ 
ent  temperatures.  The  onset  of  non-linear  behavior  occurs  at 
higher  strain  rates  for  higher  temperatures  and  lower  viscosi¬ 
ties.  The  slow  increase  in  viscosity  to  its  Newtonian  value, 
n  Q ,  at  the  low  strain  rates  corresponds  to  shear  relaxation 
effects  in  the  glass,  and  the  time  to  reach  the  steady-state, 
Newtonian  behavior  is  proportional  to  the  average  shear  relaxa¬ 
tion  time,  x q  *  n0/G  where  G  is  the  instantaneous  shear  modulus. 

• 

As  the  strain  rate  is  increased,  three  significant  effects 
can  be  observed.  First,  the  apparent  viscosity  approaches  a 
steady-state  value  for  long  times.  This  steady-state,  apparent 
viscosity  decreases  with  increasing  strain  rate.  Second,  there 
is  a  short-time  overshoot  of  this  steady-state  value  during  the 
shear  relaxation  period.  Third,  the  slower  shear  relaxation 
processes  appear  to  be  rhort-circuited  by  a  faster  process  which 
leads  to  the  lowered  viscosity.  These  results  are  qualitatively 
identical  to  the  behavior  of  Lennard-Jones  spheres  in  the  MD 
experiment  (Fig.  1).  This  process  appears  to  be  a  structural 
breakdown  mechanism  and  occurs  at  earlier  times  for  increased 
strain  rates  (i.e.  compare  the  curves  with  l  -  72*  ’O'Vsec.  and 
18*  10”4/sec.  in  Fig.  5).  It  appears  that  the  high  strain  rates 
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cause  a  breakdown  of  the  glass  structure  which  in  turn  allows 
stress  relaxation  to  occur  by  the  faster  available  mechanisms 
rather  than  the  sum  of  all  relaxation  mechanisms.  Because,  the 
Newtonian  shear  relaxation  processes  and  the  structural  break¬ 
down  processes  have  different  relaxation  times  and  amplitudes, 
there  is  an  overshoot  in  the  apparent  viscosity. 

Bach  set  of  tests  shown  on  Figs.  5  and  6  was  conducted  on  a 
single  fiber.  The  measurements  gave  the  same  steady-state  vis¬ 
cosity  values  at  each  given  2  whether  the  measurements  were  made 
following  a  higher  or  a  lower  i.  Therefore,  the  steady-state 
viscosity  values  appear  to  be  independent  of  past  history.  In 
each  figure,  the  highest  c  value  corresponds  to  &  failure  of  the 
fiber.  In  Fig.  5,  the  viscosity  appears  to  be  approaching 

steady-state  when  failure  occurs. 

* 

B.  Heating  Effects 

It  is  apparent  that  the  observed  non-linear  behavior  of  the 
silicate  glass  fibers  is  similar  to  the  behavior  of  the  Lennard- 
Jones  glass.  Before  a  discussion  of  similarities  ard  differ¬ 
ences  between  these  two  results,  it  is  necessary  to  establish 
that  test  conditions  were  similar.  For  example,  in  the  MD 
experiments,  the  temperature  was  maintained  constant.  Therefore 
it  is  necessary  to  determine  whether  the  silicate  glass  fibers 
were  heated  above  the  furnace  temperature  at  the  higher  strain 
rates,  and  whether  the  heating,  if  present,  would  be  sufficient 
to  cause  the  decrease  that  was  observed  in  the  viscosity.  For 
example,  the  drop  in  observed  viscosity  at  563°C  from  the  Newto¬ 
nian  value  of  3.2  *  1012P  to  the  pseudo-plastic  value  of  1.1  x 
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1012P  (at  the  strain  rate  of  7.2  *  lO^sec”1)  corresponds  to  a 
temperature  increase  of  10*C. 

A  lower  limit  to  the  rate  of  heat  loss  from  a  fiber  can  be 
calculated  using  the  radiative  heat  transfer  coefficient,  H, 
defined  for  small  temperature  differences  between  the  fiber  and 
the  furnace  ast 

H  «  4eST3  (3) 

where  e  is  the  glass  emissivity,  S  is  the  Stefan-Boltzmann  con¬ 
stant  and  T  is  the  absolute  temperature.  The  rate  of  heat  loss 
dQ/dt  is  then  given  by 

dQ 

^-HAsAT  (4) 

where  As  is  the  surface  area  of  the  fiber  in  the  furnace  and  AT 
is  the  temperature  difference  between  the  fiber  and  the  furnace. 
Using  S  «  5.67  x  10“® J/sec*m2  •  & ,  T  ■  836K#  and  e  ■  0.9  we  cal¬ 
culated  a  value  for  H  ■  1.2  x  102J/K«m2*sec.  This  value  is  an 
underestimate  of  the  total  heat  transfer  coefficient  since  only 
radiative  mechanisms  have  been  considered.  Paek  and  Kurkjian® 
have  estimated  H  from  cooling  rate  measurements  made  on  glass 
fibers  with  similar  dimensions  as  those  tested  here.  Their 
estimate  yielded  H  ■  2.9  x  102 J/K»m2 *sec  which  is  consistent 
with  our  calculation. 

A  calculation  of  the  rate  of  doing  work  on  the  fiber  was 
made  to  determine  an  upper  limit  on  any  resulting  temperature 
rise.  The  rate  of  doing  work  is 


dW 

dt 
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=  a  e  V 


where  V  is  the  volume  of  the  fiber  given  by  4*A(t).  Under 
steady-state  conditions  the  maximum  rate  of  heat  input  due  to 
conversion  of  all  the  work  to  heat  was  calculated  to  be  3.37  x 
10“4J/sec.  Since  the  surface  area  of  the  fiber  in  the  furnace 
was  approximately  1cm2  the  largest  temperature  rise  that  could 
be  supported  considering  only  radiative  heat  loss  is  less  than 
0.03°C.  This  is  far  below  the  value  of  10°C  necessary  to  cause 
the  observed  non-linear  behavior.  Therefore,  we  can  reasonably 
conclude  that  the  non-linear  behavior  of  the  glass  during  these 
experiments  was  not  a  result  of  localized  heating  of  the  fiber, 
and  the  test  does  duplicate  the  isothermal  conditions  of  the  MD 
experiments. 

C.  Mechanism 

During  the  MD  experiments,  it  is  possible  to  arrest  the 
system  and  examine  its  structure.  The  examinations  of  systems 
having  undergone  extensive  non-linear  behavior  showed  definite 
structural  changes  indicative  of  a  layering  effect  for  planar 
shear.  The  layers  appear  to  form  almost  parallel  to  the  shear 
planes  with  a  small  angular  deviation  from  the  shear  direction. 
The  angles  decrease  with  increasing  shear  rates.2  Real  glasses 
may  also  be  suddenly  frozen  by  rapid  cooling  to  temperatures 
below  the  glass  transition  temperature.  However,  structural 
examinations  of  real  glasses  on  the  molecular  level  are  diffi¬ 
cult  and  cannot  be  made  directly.  Therefore,  we  have  not  yet 
attempted  to  gather  structural  correspondence  tor  the  non-linear 
behavior. 
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lt  is  possible  however/  to  discuss  the  non-linear  behavior 
of  real/  inorganic  glasses  by  comparison  with  MD  results.  The 
reduced  steady-state  viscosity  (apparent  viscosity  divided  by 
Newtonian  viscosity,  n/n0  )  can  be  plotted  versus  stress  or 
versus  reduced  strain  rate  (strain  rate  times  the  Newtonian 
average  relaxation  time,  et^)  for  each  temperature.  In  both 
instances,  the  data  reduces  to  the  same  general  shape,  although 
it  appears  to  us  that  the  reduction  with  reduced  strain  rate  is 
better.  Both  inorganic  glass  experiments  and  the  MD  data  follow 
the  same  general  behavior. 

The  decrease  in  viscosity  with  increasing  strain  rate 

results  from  an  asymptotic  approach  of  the  sustained  steady- 

state  stress  to  a  maximum  value  (<»iimit^  at  very  high  strain 

rates.  The  existence  of  a  limit  in  the  sustainable  steady-state 

• 

stress  indicates  that  if  the  system  is  placed  under  a  stress 
greater  than  the  limit,  steady-state  conditions  cannot  be  main¬ 
tained  and  catastrophic  failure  ensues.  This  stress  limit, 
therefore,  can  be  interpreted  as  the  actual  cohesive  strength  of 
the  material.  This  result  shows  a  unique  and  valuable  feature 
of  the  non-Newtonian  viscosity  studies  since  the  actual  cohesive 
strength  of  the  material  is  obtained  here  without  a  need  to 
fracture  the  material. 

A  calculation  of  the  limiting  stress  requires  some  extrapo¬ 
lation  of  the  data  to  very  high  shear  rates.  In  this  task,  we 
rely  on  the  similarity  between  the  measurements  conducted  on  the 
soda-lime  silica  glass  and  the  MD  calculations  on  the  Lennard- 
Jones  glass.  The  normalized  viscosity  of  both  systems  reduces 
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to  the  same  curve  when  the  shear  rate  is  normalized  through  the 
function  erQ/a  where  o  is  an  adjustable  parameter.  The  value  of 
o  is  found  to  be  0.0063  for  the  soda-lime  silica  glass,  0.0053 
for  the  rubidium  silicate  glass  and  0.07/  for  the  Lennard-Jones 
glass.  Reduced  in  this  fashion,  all  data  follow  the  same  beha¬ 
vior,  independent  of  temperature  and  composition  as  shown  in 
Fig.  7. 

The  fall-off  in  reduced  viscosity  at  increasing  strain 
rates,  thus  appears  to  result  from  dynamic  changes  in  the  struc¬ 
ture  of  the  material  which  accompany  the  large  applied  strain 
rates.  Reduction  of  this  data  from  different  temperatures  on 
the  same  glass  by  use  of  the  average  Newtonian  shear  relaxation 
time  indicates  that  the  structural  rearrangement  is  controlled 
by  the  shear  flow  processes  in  the  glass. 

The  functional  dependence  of  the  decrease  in  normalized 
viscosity  can  be  obtained  from  an  equation  based  on  the  concept 
of  a  limiting  stress  used  by  Bair  and  Winer  to  discuss  similar 
behavior  for  highly  viscous  organic  lubricants.9  The  equation 
was  derived  by  Montrose  using  semi-empirical  arguments. The 
equation  yields  a  simple  relationship  between  reduced  viscosity 
and  the  normalized  strain  rate  function,  eij/oi 


n0  1  +  e i Q/a 

where  the  factor  a  is  found  to  be  the  ratio  of  the  maximum 
stress  sustained  by  the  system  and  the  instantaneous  shear 
modulus  of  the  glass,  The  fit  of  this  equation  to  the 

data  is  also  shown  in  Fig.  7. 


-1.3- 


0 1  im i t  tfie  maxi,num  stress  developed  in  the  system  under 
steady-state  conditions  as  the  strain  rate  goes  to  infinity 
without  material  failure.  Therefore,  0  limit  represents  the 
actual  cohesive  strength  of  the  material.  This  actual  cohesive 
strength  is  interpreted  as  the  maximum  stress  sustainable  by  the 
material  when  the  stress  is  applied  at  a  rate  slower  than  the 
effective  relaxation  time  of  the  glass.  Its  existence  suggests 
that  when  viscous  or  plastic  flow  occurs  under  an  increasing 
tensile  or  shear  stress,  there  is  a  stress  value  where  steady- 
state  flow  cannot  be  sustained  and  the  material  fractures  at  the 
point  where  the  stress  is  applied.  In  the  case  of  loaded 
fibers,  this  result  suggests  that  as  the  diameter  of  the  fiber 
undergoes  localized  necking,  a  point  is  reached  when  oappiiecj  > 
°limit  and  fracfcure  occurs,  rather  than  continued  necking  of  the 
fiber. 

In  the  MD  experiments,  an  independent  measurement  of 
°limit/G  was  conducted  and  agrees  very  well  with  the  value  of  o 
obtained  by  fitting  the  non-linear  viscosity  to  Eq.  6. 

This  generalized  behavior  of  inorganic  glasses  at  different 
temperatures  and  composition  is  an  encouraging  basis  for  the 
development  of  molecular  models  to  interpret  non-linear  behavior 
in  inorganic  glasses.  The  strong  similarities  to  the  MD  glasses 
indicate  that  the  effect  is  indeed  a  result  of  structural  rear¬ 
rangements  in  the  glasses.  The  parameter,  a  appears  to  be 
linked  to  the  strength  of  some  average  structural  bond  of  the 
material  at  the  temperature  of  measurement.  Therefore,  this 
exp  riment  offers  a  direct  measurement  of  the  cohesive  strength 
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of  different  glasses  and  a  possibility  for  measuring  its 
temperature  dependence.  The  values  obtained  here  show  that  both 
oxide  glasses  can  only  sustain  stresses  about  1/160  of  their 
shear  modulus  while  the  Lennard-Jones  glasses  can  go  up  to  about 
1/10  of  the  shear  modulus.  This  result  is  reasonable  from 
structural  considerations  since  the  bonding  forces  in  the  LJ 
glass  are  spherically  symmetric,  while  a  significant  portion  of 
the  bonding  of  the  oxide  glasses  is  covalent  and  therefore  high¬ 
ly  directional.  It  is  expected  that  these  covalent  bonds  are 
much  less  resistant  to  shear  stresses  since  they  can  break  by 
bond  rotation  than  are  the  non-directional  ionic  bonds  which 
break  by  bond  extension  alone.  Estimating  a  shear  modulus  for 
the  soda-lime  silica  glass  of  2.2  x  10loPa,  based  on  high  tem¬ 
perature  ultrasonic  sound  velocity  measurements,  our  interprets- 

% 

tion  of  the  reported  measurements  yield  an  actual  cohesive 
strength  of  1.4  x  108Pa  (20  Kpsi)  for  this  glass  at  temperatures 
near  560#C,  and  Newtonian  viscosities  of  101 °-101 xPa» s  (1011- 
101 2Poise) . 
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SUMMARY  AND  CONCLUSIONS 

The  viscosity  of  a  stable,  soda-lime  silica  glass  was  meas¬ 
ured  at  high  shear  rates.  The  results  showed  a  non-Newtonian 
viscosity  behavior  of  the  pseudo-plastic  type.  As  a  function  of 
time,  the  viscosity  at  high,  constant  strain  rates  first 
increased  to  a  maximum,  then  decreased  to  a  constant  time-inde¬ 
pendent  value,  lower  than  the  expected  Newtonian  viscosity.  The 
time  independent  behavior  at  long  times  yields  a  steady-state 
viscosity  which  is  a  function  of  the  applied  sttain  rate. 

The  measurements  when  compared  to  molecular  dynamic  calcu¬ 
lations  on  a  Lennard-Jones  material  show  unexpected  similarity 
in  the  time  and  temperature  dependence  of  the  non-Newtonian  vis¬ 
cosity.  The  observed  behavior  also  agrees  with  earlier  work  on 
a  rubidium  silicate  glass.  Comparison  of  these  data  indicates 
that  the  time  evolution  of  the  viscosity  of  the  silicate  glasses 
corresponds  to  that  calculated  for  the  Lennard-Jones  glasses 
which  is  known  to  result  from  a  structural  rearrangement  in  the 
material  under  high  deformation  rates. 

.  An  examination  of  the  combined  data  of  the  inorganic 
glasses  and  the  ND  glass  showed  that  the  steady-state,  sustained 
stress  under  applied  constant  strain  rates  approaches  a  limiting 
stress  value,  for  infinite  strain  rates.  This  limiting 
sustained  stress  is  interpreted  as  the  actual  cohesive  strength 
of  the  glass  since  any  applied  stress  greater  than  cannot 
be  sustained  under  steady-state  conditions  and  leads  to  catas¬ 
trophic  failure.  The  extrapolation  of  data  using  a  semi-empiri- 
cal  model  showed  that  0  limit  for  the  two  inorganic  glasses 
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cinalyzed  is  near  )/l60  times  the  rigid  shear  modulus,  or 
approximately  20Kpsi  for  the  soda  lime  silica  glass.  These 
results  suggest  a  mechanism  for  material  failure  at  high  temper¬ 
atures  (under  conditions  of  plastic  flow).  By  this  mechanism, 
an  applied  shear  stress  greater  than  o^m^t  induces  initially  a 
plastic  flow  reaction  which  changes  to  fracture  when  the  struc¬ 
ture  attempts  to  adjust  to  the  stressed  condition.  An  applied 
shear  stress  lower  than  induces  a  plastic  flow  reaction 
which  reaches  a  steady  state  strain  rate  as  the  structure  rear¬ 
ranges  itself  to  accomodate  the  stress. 

ACKNOWLEDGMENT 

The  molecular  dynamics  experiments  described  in  this  paper 
were  supported  by  the  Office  of  Naval  Research  (Contract  No. 
N00014-75-C-0856 ) . 


120 


-17- 


FIGURE  CAPTIONS 


Fig.  1  --  Results  from  molecular  dynamics  calculations.  Shear 
stress  versus  time  for  the  shear  strain  rates  shown  at 
the  right,  switched  on  at  time  zero.  Note  the  stress 
overshoot  at  the  upper  strain  rate.  (See  Ref.  2  for 
details) . 

Fig,  2  —  Details  of  the  furnace  and  load  cell  assemblies.  The 

end  heaters,  1  and  3  were  adjusted  to  obtain  a  uniform 
temperature  in  the  furnace. 

Fig.  3  —  Variation  of  stress  versus  strain  rate  for  the  soda- 
lime-silica  glass.  The  deviation  from  Newtonian  beha¬ 
vior  is  clearly  seen  for  strain  rates  above  1.5  x 

10“4sec~1  . 

Fig.  4  —  Newtonian  Viscosity  of  the  measured  soda-lime-silica 
glass. 

Fig.  5  —  Measured  viscosity  versus  time  for  the  strain  rates 

shown,  switched  on  at  time  zero.  The  dashed  line  is 
the  Newtonian  value. 

Fig.  6  —  Measured  viscosity  versus  time  for  the  strain  rates 

shown,  switched  on  at  time  zero.  The  dashed  line  is 
the  Newtonian  value. 

Fig.  7  —  Plots  of  viscosity  reduced  by  Jthe  Newtonian  value 
versus  normalized  strain  rates,  et.  for  various  tem¬ 
peratures.  The  solid  line  represents  equation  (6)  and 
shows  how  a  reduction  of  the  strain  rate  by  the  single 
parameter  of:  G/o  (limit)  can  fit  the  data  from  three 
totally  different  materials: 

a)  Rubidium-silicate  data  from  Li  and  Uhlmann4  ana¬ 
lyzed  by  the  method  presented  in  this  paper.  The  tem¬ 
peratures  are  as  follows:  open  circles  »  528#C,  solid 
trianges  »  555#C,  open  triangles  ■  536#C,  squares  = 
501*0 ,  solid  circles  ■  480°C. 

b)  Soda-lime-silica  glass  whose  measurement  is  des¬ 
cribed  here.  The  temperatures  are  as  follows:  cir¬ 
cles  ■  563°C,  diamonds  ■  574#C,  triangles  *  593°C,  and 
squares  ■  596°C.  The  two  solid  points  are  estimates 
of  viscosity  before  failure. 

c)  Molecular  dynamics  calculations  on  a  Lennard-Jones 

glass.2  Here  different  densities  are  used  to  repre¬ 
sent  different  thermodynamic  states.  The  densities 
are  in  terms  of  the  triple  point  density,  p™:  circles 
*  lpT,  triangles  =■  l.lpT,  diamonds  »  1.2pT,  and 

squares  =  1.5pT. 
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